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The Inverse Energy Cascade of Two-Dimensional Turbulence 

Michael K. Rivera, Ph.D. 
University of Pittsburgh, 2000 

This thesis presents an experimental study of the inverse energy cascade as it occurs 
in an electromagnetically forced soap film. It focuses on characterizing important 
features of the inverse cascade such as it's range, how energy is distributed over the 
range and how energy flows through the range. The thesis also probes the assumption 
of scale invariance that is associated with the existence of an inverse cascade. These 
investigations demonstrate that the extent of the inverse cascade range and the be- 
havior of the energy distribution are in agreement with dimensional predictions. The 
energy flow in the inverse cascade range is shown to be well described by exact math- 
ematical predictions obtained from the Navier-Stokes equation. At no time does the 
energy flow in the inverse cascade range produced by the e-m cell behave inertially or 
in a scale invariant manner. Evidence that the cascade could become scale invariant 
should an inertial range develop is presented, as are the requirements that a system 
must satisfy to create such an inertial range. 
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Chapter 1 
Introduction 



One of the curious predictions in turbulence theory is that there might possibly exist 
a range of length scales in a two-dimensional (2D) turbulent fluid over which kinetic 
energy is transferred from small to large length scales. That this range could exist was 
first predicted in the late 1960's by Kraichnan[l]. Numerical simulations that followed 
yielded varying degrees of agreement with this prediction[2, 3, 4]. Experimental 
verification of the existence of such a range did not come about until much later due 
to the difficulty inherent in building and maintaining a system which approximates a 
2D fluid[5, 6, 7]. This thesis presents an investigation of the 2D inverse energy cascade 
in a new apparatus, the electromagnetically forced soap film. What follows in this 
chapter is a description of the phenomenology surrounding the inverse energy cascade 
and a discussion of the experiments that have attempted to probe it's properties. 



1.1 The Inverse Energy Cascade 

The phenomenology of turbulence, in three-dimensions (3D) or 2D, is usually phrased 
in terms of "eddies". An eddy itself is not a well defined object, though there have 
been many recent attempts using wavelets to better define the concept [8]. Loosely 
speaking it is a region in a fiuid that is behaving coherently. The extent of an eddy is 
dictated by boundaries within which an arbitrary determination is made that some 
sort of structure exists. Thus an eddy can be a single large region of rotation, such 
as the whirlpool which forms above a bathroom drain. Or an eddy can be a large 
region containing many smaller eddies which are interacting with one another while 
behaving distinctly (again by an arbitrary determination) from other neighboring 
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Figure 1.1: Two pictures of the "eddy" concept for a 2D fluid: (a) a single large eddy 
and (b) a large eddy made from many interacting smaller eddies. 

clusters of eddies. These two ideas are drawn in Fig. 1.1 for the case of a 2D fluid. 
Eddies in 3D are much more difficult to picture. Two important properties that are 
associated with an eddy are size and energy. These two properties allow predictions 
about energy motion in fluids to be made if some knowledge of how eddies interact 
in the system is known. 

In 2D fluids, one way in which eddies (which assume the more familiar label 
"vortices" in 2D) interact with each other is through a process known as "vortex 
cannibalization" . A cannibalization event is when two neighboring eddies of like 
rotational sense merge to form a single larger eddy. When cannibalization occurs 
energy flows out of the length scales of the initial eddies and into the length scale 
of the flnal eddy. Since the flnal eddy is larger than the initial ones, cannibalization 
results in the flow of energy from small to large length scales. 

In a 2D turbulent fluid, many eddies are generally created at a small length scale 
called the energy injection scale, Vinj. The expectation is that through interaction 
by cannibalization these small eddies cluster and merge into larger eddies. These 
larger eddies are also expected to cluster and merge to form even larger eddies and so 
on. This means that energy, initially injected into the turbulence at the length scale 
Tinj should gradually be moved by consecutive cannibalization events to larger length 
scales. This type of energy motion constitutes an inverse energy cascade[l]. 

Using the eddy concept has the advantage of highlighting two important features 
associated with the existence of an inverse energy cascade: scale invariance and local- 
ity of interaction. The flrst of these can be understood by looking again at Fig. 1.1(b) 
which shows many smaller like signed eddies clustering to form a single larger eddy. 
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Presumably, the small eddies in the figure are themselves formed by the clustering 
of even smaller eddies, which in turn are formed by even smaller eddies. Likewise, 
the large eddy cluster in the figure is most likely interacting with other eddy clusters 
in the system. As long as the eddies at the very smallest scale, the injection scale, 
are being continuously created to replenish those which are cannibalized, the inverse 
cascade range is scale-invariant. That is to say that no length scale in the inverse 
cascade range can be distinguished from any other length scale that is also in the 
range. Scale invariance is exceedingly important from a theoretical stand point. The 
assumption of scale invariance of fields, such as the probability of velocity difference 
on a length scale r, allow important predictions about turbulence to be made (see 
chapter 5) [9]. 

Before discussing locality, a delicate point must be made. If the eddies at the 
injection length scale are not being continuously replenished then the number of 
eddies at the smallest scales gradually begins to decrease as more and more eddies 
are lost to cannibalization events. To maintain an inverse cascade range, then, the 
turbulence has to be continuously forced. That is, eddies must be continuously created 
at the energy injection scale. If the turbulence is not forced then the cascade range 
will eventually consume itself from small scales up, ultimately leading to a state 
which can be described as a diffuse gas of large individual eddies (eddies not made of 
clusters of smaller eddies) [10, 11]. The term "coarsening" is used to describe decaying 
2D turbulence's behavior in order to distinguish it from the inverse energy cascade. 

The second property assumed to hold in the inverse cascade is locality of interac- 
tion. This property refers to constraints on the manner in which eddies interact. If 
an eddy of very small size is close to, or embedded in, an eddy of exceedingly large 
size, the small eddy will merely be swept along by the large eddy and not strongly 
deformed. Likewise the large eddy will not be significantly effected by it's small com- 
panion. Since neither of the eddies is strongly deformed, the cannibalization process 
is expected to happen over a long period of time, if at all[9]. On the other hand, two 
neighboring eddies of similar size interact and deform one another strongly, and thus 
the cannibalization happens swiftly. Energy transfer by cannibalization is therefore 
most efficient when occurring between similarly sized length scales; this is what is 
meant by locality. Due to locality, the kinetic energy at small scales in the inverse 
cascade is expected to be moved to large scales in a continuous manner, stepping 
through the intervening length scales by local interactions rather than making large 
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length scale jumps by the merger of a small and large eddy. Hence the term cascade. 

The picture of 2D turbulence and it's inverse cascade is now almost complete. 
Energy is continually injected into a fluid in the form of small eddies. These small 
eddies cluster to form large eddies moving energy to larger scales. In turn the eddy 
clusters themselves cluster to form larger clusters of eddy clusters, etcetera. It is 
the etcetera that is of concern. At what point does the vortex merger process and 
growth of larger and larger eddies stop? That is, how is the energy injected into a 
2D turbulent fluid thermostated? 

Consider first the thermostating mechanism in 3D turbulence. In 3D turbulence 
there exists a direct energy cascade, where instead of energy being moved from small 
to large scales by eddy merger the opposite happens; energy is moved from large 
to small scales by eddy stretching (commonly called vortex stretching). Eventually, 
through continuous vortex stretching, a smallest eddy scale is reached, at which point 
the kinetic energy contained in these small eddies is dissipated into heat by the fluids 
internal viscosity. All of the energy that is injected into the large length scales of a 
3D turbulent fluid is eventually exhausted by viscosity at small length scales[9]. 

Internal viscosity is a short range force, only becoming a good thermostat when 
the kinetic energy reaches small length scales[9]. Thermostating is not an issue in 
3D where the direct cascade takes energy down to such small scales. In 2D, however, 
the inverse cascade moves energy away from small scales. Therefore viscosity has 
no chance to exhaust the injected energy. An ideal 2D turbulent fluid driven to a 
state of turbulence with a continuous forcing would never be in a steady state since 
the total energy in the flow would continue to build up as larger and larger eddies 
form [12]. What is needed to maintain 2D forced turbulence in a steady state and stop 
the inverse cascade process is some sort of external dissipation mechanism which is 
an effective thermostat at large length scales. In other words, some sort of dissipation 
mechanism that is not internal to the fluid itself must exist to take energy out of large 
length scales and dictate the largest size eddies that can be formed by the cascade. 

Fortunately, 2D experiments are almost always coupled to the surrounding 3D 
environment by frictional forces[13, 14, 7]. In these experiments this external fric- 
tional force provide the turbulence with an effective large scale thermostat and sets 
the largest length scales which can be reached by the inverse cascade process. The 
inclusion of an external thermostat completes this phenomenological description of 
the inverse energy cascade. 
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1.2 History of Experiments 

Laboratory experiments which have attempted to probe the inverse energy cascade 
of 2D turbulence fall into two major categories: soap films and stratified shallow 
layers of fluid. The soap film experiments in 2D fluid mechanics were initiated by 
Couder in the early 1980's[15]. This early work investigated coarse features of both 
2D turbulence and 2D hydrodynamics. Further attempts at using the soap film to 
measure 2D turbulence in the search for an inverse cascade were done by Gharib and 
Derango a few years later [16]. The experimental system that was used by Gharib 
and Derango, called a soap film tunnel, was later perfected by Kellay et a/. [17] and 
Rutgers et a/. [14]. 

The soap film tunnel is a 2D equivalent of the wind tunnel, which is the mainstay 
of 3D turbulence research. The manner in which turbulence is created in each is 
identical. A grid, or some other obstacle, is placed in the path of a swiftly moving 
mean flow. If the flow speed is fast enough, the fluid becomes turbulent downstream 
from the grid. Such 2D and 3D tunnels are shown in Fig. 1.2. Soap film tunnels would 
seem to be ideal 2D fluids for performing turbulence research in because their aspect 
ratios are exceedingly large (many cm across to a few micrometers thick) and thus the 
fluid flow is almost entirely two-dimensional. There are, however, difficulties inherent 
in the use of soap films. For example thin films couple strongly to the air and the 
magnitude of their internal viscosity is large compared to that of water. For the most 
part these difficulties are thought to be mitigated by clever experimental techniques, 
such as the use of vacuum chambers[14] or by using thick lO/im) films[ll]. 

Every attempt to study the inverse energy cascade in a 2D soap film tunnel with 
the configuration shown in Fig. 1.2 has met with failure. This is not a disparaging 
comment about the researchers involved in the effort. Indeed their considerable skill 
eventually tamed the delicate and whimsical soap films into a useful experimental 
system. The lack of inverse cascade in these systems reflects the fact that the con- 
figuration shown in Fig. 1.2 creates decaying 2D turbulence. The eddies that are 
injected at the grid are not replenished as the fluid moves downstream. By the dis- 
cussion in the last section this means that the system does not have the ability to 
form the eddy clusters that is expected of an inverse cascade. 

Once the understanding that 2D turbulence needs to be forced for an inverse 
cascade to be present was reached, the fllm tunnel design was modifled to create 
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Figure 1.2: A 3D wind tunnel creating turbulence and it's 2D equivalent: the soap 
film tunnel 

forced turbulence [6]. The design of the film tunnel is identical to that shown before 
except for the orientation of the turbulence producing grid. Instead of having a single 
grid oriented perpendicular to the flow direction, two grids were oriented at angles 
to the flow direction so that they formed two sides of a triangle with the tip of the 
triangle oriented upstream. This modified form is shown in Fig. 1.3. An area of 
forced turbulence exists in the interior between the two grids since it is here that 
vortices created at the grids are able to diffuse into the interior and replenish those 
lost to cannibalization. 

Though certain properties of the inverse cascade can be investigated with such 
a modified film tunnel, the setup is not ideal because the turbulence it creates is 
inhomogeneous. One can imagine that the fiuctuation in the driven turbulence area 
near the grids are quite large compared to those in the interior. Homogeneity is a 
critical simplifying assumption in almost all areas of turbulence theory. The inhomo- 
geneity of the modified film tunnel, then, has devastating consequences with regard 
to comparing results with theory. A more ideal setup would involve the injection of 
vortices directly into the body of the fiuid by some sort of external force, rather than 
injecting from the boundaries. This is a difficult task to do in soap film tunnels. 

A system which does achieve such an injection of vortices falls into the second class 
of 2D turbulence experiments: stratified shallow layers of fiuid. For the most part, the 



Chapter 1. Introduction 



7 



Figure 1.3: A 2D soap film tunnel creating forced 2D turbulence. 

use of such layers in turbulence research has been pioneered by Tabeling and others 
in the early 1990 's. A successful observation of an inverse cascade regime was reached 
a few years later by the same group [18]. The stratified shallow layer apparatus, as 
it's name suggest, suspends a layer of pure water, above a higher density layer of salt 
water. The layers in question are only a few millimeters thick, and the area tends to 
be on the order of ten centimeters so that the aspect ratio, while not nearing that of 
soap films, is still large. The salt water layer is subject to a current flowing in it's 
plane, and placed in a spatially varying magnetic field. The resultant Lorentz force 
acts directly on the fluid layer driving it to a state of turbulence. Stratification helps 
to impose a measure of two- dimensionality to the fluctuations, thus one has 2D forced 
turbulence. 

Note that in this system the force is acting directly on the fluid, unlike the modi- 
fied film tunnel where forcing happened near the grid. This restores homogeneity to 
the system, allowing accurate comparison with theory. The place that shallow layers 
suffer is in their approximation of a two-dimensional fluid. The bottom of the fluid 
container in shallow layer systems enforces a no slip boundary on the lower surface 
of the fluid. If the velocity in the fluid becomes large, a strong shearing can develop 
between the upper and lower layers of fluid. This inevitably causes mixing which 
destroys the stratification and sacrifices two- dimensionality. Thus, shallow layer sys- 
tems are severely limited in the strength of the turbulent fluctuations which they can 
successfully explore. 
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This brief experimental history, then, shows that there have been two systems used 
to investigate the inverse cascade, neither of which are ideaL However, each system 
complements the others difficulties almost perfectly. Film tunnels are great 2D fluids, 
but not easily forced homogeneously. Stratified layers are marginal 2D fiuids, but 
easily forced homogeneously. What is needed then, as an ideal test apparatus for the 
inverse cascade, is a combination of these two systems that takes advantage of their 
benefits without inheriting their difficulties. What is needed is the electromagnetically 
forced soap film, simply called the e-m cell. 

1.3 Thesis Overview 

This thesis presents an experimental study of the inverse energy cascade as it occurs 
in an electromagnetically forced soap film. In particular it focuses on characterizing 
important features of the inverse cascade such as it's range, how energy is distributed 
over the range and how energy fiows through the range. The thesis also probes the 
assumption of scale invariance that is associated with the existence of an inverse 
cascade. 

Chapter 2 describes the workings of the e-m cell and measurement apparatus. The 
basic design and implementation of the e-m cell is reviewed in the first two sections. 
How the frictional coupling of the soap film to the air is controlled is described in the 
following section. The fourth section explains certain limitations on apparatus size 
that are imposed by the existence of gravity. This size constraint is one of the chief 
difficulties that limit results throughout the rest of the thesis. The fifth section in 
the chapter is an overview of the measurement system that was developed to extract 
velocity information from the e-m cell. The chapter is concluded with a brief overview 
of the operation of the e-m cell. 

Chapter 3 is a systematic check that the e-m cell does behave as a 2D fiuid. 
The first section motivates the need for this test and the second derives the relevant 
mathematical relationships necessary for such a test. The third section compares the 
data from the cell to these derived relationships to verify that the e-m cell behaves as a 
2D fiuid. This section also establishes a model for the external dissipation mechanism 
in the e-m cell. A self consistency check is also performed to help strengthen this 
verification. 

Chapter 4 begins the systematic investigation of the properties of the inverse 
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energy cascade. The extent of the range and how energy is distributed over this 
range is measured and presented in the first section. The manner of energy transfer 
is described in the following section. In particular, this second section attempts to 
determine if the energy flow is "inertial" , a property which is necessary if one wants 
results which are universal to all 2D turbulence systems. 

Finally, chapter 5 attempts to determine if the inverse cascade is scale invariant. 
Predictions for structure functions in a scale invariant fluid are presented in the first 
section. The next two sections attempt to extend results from the e-m cell to test these 
predictions, as well as explain why one should be critical of such extensions. While 
conclusions presented in earlier chapters are quite strong, the conclusions drawn in 
this chapter are weak. This weakness stems from the size limitations presented in 
chapter 2 and can not be overcome in the current apparatus. 



Chapter 2 



The E-M Cell 

A fluid which carries a current density, J, in the presence of a magnetic field, B, is 
subject to a force per unit mass F = J x B which drives fluid motion. This principle 
has been used in earlier experiments to excite motion in shallow layers of electrolytic 
fluid[5, 19], and the techniques used in these experiments may be readily adapted 
for use in soap films. The result of such adaptation is the electromagnetically forced 
soap film, called briefly an e-m cell. The e-m cell is a useful tool in the study of 2D 
hydrodynamics, and in particular 2D turbulence [7]. This chapter reviews the design 
and operation of the e-m cell, as well as the measurement technique used to obtain 
velocity data from it. 

2.1 The E-M Cell 

The main component of the e-m cell is a free standing soap film drawn from an 
electrolytic soap solution across a square frame. The frame has two opposing sides 
made from stainless steel, while the remaining sides are constructed from plastic or 
glass. A voltage difference applied to the stainless steel sides results in a current 
which lies in the plane of the film. A spatially varying external magnetic field is 
then created and oriented so that it penetrates the film plane perpendicularly. The 
resulting Lorentz force lies in the plane of the film and drives the fluid motion. A 
diagram of this is shown in Fig. 2.1. 

The electrolytic soap solution which the film is drawn from is made from 400 ml 
distilled water, 80 g ammonium chloride salt, 40 ml glycerol and 5 ml commercial 
liquid detergent (regular Dawn or Joy). To this solution particles are added up to a 
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Figure 2.1: Basic operation of the e-m cell. 



volume fraction of about 10""^. These particles, either 10 yum hollow glass spheres or 
lycopodium mushroom spores which appear as ~ 40 /im particles, are of a density 
comparable to the soap solution so that they closely follow the surrounding flow. 
These particles will be used in the measurement technique to be described later. The 
salt used to make the solution electrolytic, ammonium chloride, was chosen after 
numerous trials which included sodium chloride and potassium chloride. Ammonium 
chloride was chosen because high concentrations could be reached with little effect 
on the stability of the film. This was found not to be the case with potassium or 
sodium chloride. High concentrations of salt are necessary to keep the films electrical 
resistance as small as possible. This limits film evaporation due to Joule heating 
caused by the driving current, a necessity if a film is to be studied for long periods of 
time. 

The soap film, once drawn, is maintained at a thickness of around 50 /im. Earlier 
soap film work tended to use thin films with thicknesses of 5 /im or less. There are a 
number of advantages that come from using thick films. First, for the e-m cell, it is 
desirable to keep the electrical resistance low, again to limit evaporation due to Joule 
heating. The larger the cross sectional area through which the current is passed, then, 
the better. Another reason to use thick films is that they lose a smaller percentage of 
their energy than thin films to frictional rubbing against the surrounding air. Recall 
from chapter 1 that 2D driven turbulence requires an external dissipation mechanism 
to maintain energy balance. Air friction plays this role in the e-m cell. If the air 
friction is strong then it can easily dissipate large amounts of energy. Therefore 
the energy injected from the electromagnetic force must be exceptionally high to 
maintain a state of strong turbulence. This would require large currents which would 
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enhance the Joule heating and evaporation of the film. Using thick films reduces 
this coupling to the air and allows strong turbulence to be maintained for reasonable 
values of the driving current. A final reason for using thick films arises from the 
observation that a soap films kinematic viscosity, is dependent on it's thickness. 
Trapeznikov predicted that the effective kinematic viscosity should depend on the 
thickness, h, a.s u = i^buik + '^'^surf/h, where Uhy^ik is the kinematic viscosity of the soap 
solution and Ugurf is the 2D viscosity due to the soap film surfaces[20]. The surface 
viscosity was recently measured^ to be around I'surf ~ 1-5 x 10~^ cm^/s for soap films 
similar to that in the e-m cell[21]. Since the soap solution's viscosity is .01 cm^/s the 
effective viscosity of a 50 fim thick soap film is approximately 0.016 cm^/s. Using 
thick films, then, reduces the amount of energy lost to viscous dissipation, and by the 
same argument given above, allows strong turbulence to be maintained for reasonable 
values of driving current. 

There are, of course, disadvantages to the use of thick soap films. One is that 
the speed of a 2D density wave in the soap film is dependent on the thickness of the 
film^. Thicker films contain more interstitial fluid than thin films, and therefore have 
a lower wave velocity due to their increased mass. Therefore, as the film becomes 
thicker it begins to be more easily compressed. That is to say that it's 2D density 
couples more strongly to the velocity field in the film. A measure of the importance 
of such compressibility effects is the Mach number M = Urms/c, where c is the density 
wave speed in the medium. In the case of a 50 fim thick film the wave velocity is 
approximately 2 m/s. Velocity fiuctuations in the e-m cell are therefore kept less 
than 20 cm/s so that M < 0.1. With such a small Mach number the system does not 
develop shock waves and behaves approximately as an incompressible fiuid. 

As noted earlier. Joule heating evaporates fiuid from the film. This causes the 
average thickness of the soap film to change over time. Since thick films were used to 
minimize this effect, thickness changes due to evaporation happen at a slow rate and 
can be balanced by injecting small amounts of fluid into the fllm. In the e-m cell, soap 
solution is injected by a syringe pump through a small needle inserted through the 
plastic side of the frame as in Fig. 2.2. It is important that the fluid be burst in, as 

^This measurement utilized a flowing soap film tunnel to analyze the shedding of vortices from 
a cylinder placed in the flow. The vortex shedding frequency is sensitive to changes in the fluids 
viscosity and through proper normalization allows I'surf to be approximated. 

2D density wave in a soap film is a thickness wave, where the film bulges or shrinks in the 
third dimension. 
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Figure 2.2: Replenishing fluid lost to evaporation. 

a brief jet, and not slowly injected. Without initial momentum, fluid builds up near 
the injection point forming a droplet which eventually pops the film. Fluid injected 
with a burst shoots to the center of the e-m cell where it is quickly mixed into the 
film replenishing the lost fluid. Though not used in the experiments reported here, 
this process can be automated by monitoring the film resistance. When the resistance 
becomes to high a small burst of fluid is shot into the film raising the thickness and 
lowering the resistance to an acceptable level. 

The square frame that the film is drawn across is limited in size to 7 x 7 cm^ in 
all of the experiments reported here for reasons which will be made clear later. Both 
the stainless steel and plastic sides are milled to a sharp edge of about 45°, though 
the stainless steel edge generally dulls over time by corrosion. Recent efforts have 
attempted to replace these edges with sheets of stainless steel and glass that are less 
than 100 fim thick. The intent of this is to limit the size of the wetting region near 
the edge of the frame, as shown in Fig. 2.3. Since the film has a finite surface tension, 
there must be a pressure jump across the film surface if it is curved. A wetting region 
induces a negative curvature near the edge of the frame, therefore the pressure inside 
the film near the frame edge is less than at the film center where there is no curvature. 
This causes fluid to be forced through the film to the edges, causing the center of 
the film to drain. Elimination of these wetting regions eliminates such drainage. 
This effect is relatively weak compared to other effects which cause drainage so using 
sheets instead of edges is a low priority, sheets being difficult to work with due to 
their delicacy. All experiments reported in this thesis use edges instead of thin sheets. 
However, the preliminary work with sheets which is not reported here has provided 
encouraging early results. 

Using a 50 /^m film made from the above ammonium chloride solution on the 7x7 
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Figure 2.3: Film curvature near an edge and a plate. 



cm^ frame results in resistance values of around 1 kQ. Driving currents of between 10 
and 45 mA are used depending on the strength of turbulence desired, strength of the 
air dissipation, and strength of externally applied magnetic field. In all turbulence 
experiments, the applied voltage oscillates at 3 Hz with a square waveform. There are 
several reasons for this, one is that the current in the film causes chlorine gas bubbles 
to accumulate on the positive electrode. By oscillating the current, the bubbles form 
at both electrodes, slowing the inevitable bubble build up which eventually invades 
the film and renders it useless. Another reason for current oscillation is to eliminate 
the formation of large vortices which form and dominate regions of the fluid containing 
a force field which is sympathetic to it's motion. Once formed, such structures induce 
a spatially varying mean flow, and thus a mean shear, rendering the turbulence in 
the e-m cell inhomogeneous.^ 

A final note: due to the configuration of the e-m cell, the current which is driven in 
the plane of the film is unidirectional, namely parallel to the plastic edges of the frame 
running from one electrode to another. For the rest of the paper the current direction 
will be assumed to lie along the y axis. By the Lorentz law, this unidirectional 

■^Oscillating the current may also eliminate any polarization of charge in the e-m cell caused by 
the net motion of ions in the fluid. At this time it is unclear whether this polarization is important 
in the e-m cell. Experiments done in both laminar and turbulent flow with D.C. forcing do not 
exhibit signs of charge polarization, such as the gradual decrease in applied force due to shielding of 
the electrodes. However, these experiments where done over short (minutes) time scale. What effect 
polarization may have over a typical experiment time scale (around fifteen minutes) is unknown 
since such efi^ects are masked by the efi^ects of evaporation. 
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Figure 2.4: Array of magnets creating the spatially varying external magnetic field. 

current creates a unidirectional electromagnetic force which lies on the x axis. Due to 
this unidirectionality, the e-m cell can not behave isotropically, that is the system is 
not rotationally invariant, regardless of the symmetry of the spatial variation of the 
magnetic field. This lack of isotropy will be an exceedingly important consideration 
later in the thesis. 

2.2 The Magnet Array 

The spatially varying magnetic field used in the creation of the Lorentz force is gen- 
erated by an array of Neodymium rare earth magnets (NdFeB) placed just below the 
film as shown in Fig. 2.4. These magnets produce quite powerful fields, typically on 
the order of 0.1 T at their surface. This field decays away from the surface, however, 
with a typical length scale of the order of the magnet size. Since the magnets are 
small, generally less than 0.5 mm, they must be brought very close to the film surface 
to generate fields strong enough to drive flow. One could use larger magnets to induce 
motion. This is undesirable, though, since the magnet size dictates the energy injec- 
tion scale, rinj- Recall from the discussion in chapter 1 that Vinj must be significantly 
smaller than the system size for an inverse cascade region to exist in an experimental 
2D turbulence system. For the e-m cell described above, the system size is the frame 
size, which is limited to 7 cm. Therefore magnet arrays with injection scales less than 
0.7 cm must be used to allow for a decade of inverse cascade range. 

Several types of magnet arrays are used in the e-m cell, each distinguished by the 
type of flow it induces at small Reynolds number (it's laminar behavior) or equiva- 
lently by the characteristics of the force field it produces in the e-m cell. The first type 
of array, and by far the most important, is the Kolmogorov array. This array is made 
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Figure 2.5: Top view of the magnet arrays which create the spatially varying external 
magnetic field in the e-m cell: (a) the Kolmogorov array, (b) the square array, (c) 
stretched hexagonal array, (d) pseudo-random array. The direction of the current J 
is shown as is the coordinate axis. 
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of alternating north-south layers of long rectangular magnets of approximately 0.3 
cm in width, as shown in Fig. 2.5 a. The magnetic field it produces in the film varies 
approximately sinusoidally in space. The magnets are oriented so that this variation 
lies in the direction of the current, causing the Lorentz force in the film to have the 
form Fx = fosin{koy) (Recall Fy = because the force is unidirectional). Two mag- 
nets make a single north-south oscillation, therefore the wavelength of the sinusoidal 
variations in the force field (which is the energy injection scale) is Vinj = 0.6 cm. The 
associated injection wavenumber is = 27r/rj„j = 10 cm~^. The laminar flow this 
forcing produces is one of alternating shear layers, a flow which Kolmogorov proposed 
investigating as an interesting model to study a fluids transition from laminar to tur- 
bulent motion (hence the name). The Kolmogorov array has several properties which 
will be of importance in later analysis. The first is that the forcing is divergence free, 
i.e. V ■ -F = 0. This property will allow pressure fields to be approximated from 
velocity fields using Fourier techniques. Another property of importance is that the 
forcing is invariant to translation along the x direction. Chapter 3 will demonstrate 
how this property may be used to obtain the energy injection rate from the forcing 
without having to explicitly measure the magnitude of the forcing /q. 

The next type of magnet array, shown in Fig. 2.5 b, is called the square array. 
It is made of round magnets oriented in a tick-tack-toe arrangement with like poled 
magnets along diagonals. There are two such arrays in the e-m cell arsenal made 
from 0.3 cm and 0.6 cm magnets. The force field created by these arrays has the form 
Fx = fo{sin{kQ{x + y)) + sin^ko^x — y))). Here the wavelength along the diagonals 
is rinj = ^/2w where w is the magnet diameter, and again ko = 2'n/rinj. The 0.3 cm 
square array has the smallest injection scale of all the arrays used with rj„j = 0.42 
cm. Unfortunately these magnets are exceptionally weak, limiting the magnitude of 
the forcing one can create with reasonable currents. Also, the forcing created by these 
arrays are not divergence free. Without the ability to account for pressure much the 
usefulness of this array is limited to testing ideas of universality (Chapter 5). 

The final two arrays are seldom used and are mentioned here only for completeness. 
The first, shown in Fig. 2.5 c, is a stretched hexagonal array, called the hex array, 
and the second, shown in Fig. 2.5 d, is a pseudo-random array. The former is 
constructed from a mixture of both 0.3 cm round and 0.6 cm round magnets. The 
result could be though of as a Kolmogorov flow with a periodic variation along every 
other magnet. The injection wavenumber associated with this array is quite large. 
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Figure 2.6: Fluid between a top plate moving with velocity U and fixed bottom plate 
produces a linear velocity profile. 

limiting it's usefulness. The laminar behavior a hex array produces is a triangular 
vortex crystal. The pseudo-random array is constructed from 0.3 cm round magnets 
placed at random on an iron sheet. The positions of the magnets where generated 
via a random number generator which attempted to maximize the magnet density. 
It was constructed in a naive attempt to obtain more homogenous turbulence in the 
e-m cell. The pseudo-random flow has no well-defined laminar flow behavior. 

2.3 External Dissipation: The Air Friction 

The air friction has already been described as the external dissipation mechanism 
in the e-m cell. Chapter 3 will demonstrate that the air friction can be adequately 
modeled as a linear friction, that is if an element of the film is moving with velocity 
u it experiences a drag force from the air of the form Fdrag = —au. The following 
discussion sets the basis for this linear drag model. 

A fluid between two parallel plates separated in the X2 direction by a small distance 
one fixed and the other moving with velocity U in the xi direction, assumes a linear 
profile of the form, ui = Ux2/d (This assumes there are no other forces acting on 
the fluid, the plates are infinite and boundary conditions are no slip at either plate). 
This is shown in Fig. 2.6. The drag force per unit area exerted on the top plate by 
the fluid is fdrag = V^\d = ^ where rj is the ordinary 3D shear viscosity of the fluid 
between the plates. In the e-m cell the role of the top plate is played by the soap 
film, while the role of the lower plate is played by the magnet array. The fluid is the 
surrounding air and the length d is just the distance of the magnets to the film. Thus 
raising the magnets closer to the film increases the strength of the air drag in the e-m 
cell. 
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The drag force must be normalized by the films 2D density so that the force per 
unit mass can be considered. The 2D density of the film is given by ph where p is 
the density of the soap solution (about that of water) and h is the films thickness. 
The force per unit mass caused by the air drag on the film is then Fdrag = 
This demonstrates a point that was made earlier in this chapter, that the frictional 
coupling to the air depends not only on the magnet-film distance but on the thickness 
of the soap film as well. Keeping the film as thick as possible, then, minimizes this 
coupling. 

There are a couple of reasons to worry about the applicability of this linear friction 
model. First is that the soap film in the e-m cell is not an infinite flat plate moving 
with a constant velocity, but has a velocity which fluctuates from point to point in 
the flow. As long as the size of these velocity fluctuations, that is the size of the 
vortices, is larger than both the magnet-film distance, d, and the film thickness, /i, 
the approximation as an infinite plate should apply. Recall from the discussion above 
that the magnetic field created by the array quickly decays within the width of a 
single magnet. Therefore the magnets are always kept within one magnet width of 
the film, otherwise the magnetic field would be too weak to drive turbulence. Since 
the magnet width also dictates the smallest vortex size, rj^j, the requirement that 
d < Tinj is always met in the e-m cell. 

The second complication is the fact that there is more than one side to the film. 
The linear drag model accounts for the drag force exerted on the lower surface of the 
soap film. The upper surface of the film is also dragging along a layer of air. The 
velocity profile of the air above the film is not a simple linear shear as it is below the 
film, but a more complex Prandtl-like boundary layer due to the absence of a second 
fixed plate. Chapter 4 will show that this causes a non-negligible correction to the 
above linear drag model. 

2.4 Gravity 

To this point there has been no discussion about the orientation of the soap film 
with respect to the earth's gravitational field. Vertical orientation, that is the film 
plane lying along the gravitational field direction, is not desirable because it would 
strongly stratify the soap film, making it thinner on top than on the bottom. This 
is tantamount to both a severe change in the 2D density and a change in the films 
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Figure 2.7: A thick film droops under the action of gravity, as shown by the dotted 
hne. A box enclosing the top of the e-m cell frame is brought to a lower pressure than 
the surrounding environment to balance gravity. 



kinematic viscosity from the top of the film to the bottom. In other words the film 
becomes an inhomogeneous fluid. Vertical orientation must therefore be discarded 
and a horizontal orientation used to allow the film to approximate a homogenous 
fluid. 

A horizontally oriented soap film droops under the force of the earth's gravitational 
field. This effect is exacerbated by the fact that the film is 50 fim thick. To balance 
gravity, a box enclosing the region on the top surface of the soap film is evacuated 
of a small amount of air, as in Fig. 2.7. This lowers the pressure in the box causing 
a pressure gradient across the film plane and thus a force opposing gravity. Enough 
pressure is drawn to almost exactly balance the gravitational field. Unfortunately the 
larger the soap film, the more delicate this balance becomes. This is the reason that 
soap films used in the e-m cell are limited to sizes under 7x7 cm. 

This pressure balance is delicate and can be disrupted by the evaporation of fluid 
from the soap film into the container above the film. It must be constantly monitored 
to make sure that the film stays at the same level. This is done in the e-m cell by 
reflecting a laser light off a portion of the film near the middle of one of the edges. 
A drooping film deflects the beam, and this deflection can be monitored by various 
techniques, for example using a position sensitive detector or a linear CCD array. A 
feedback loop based on the deflection measurement can be easily constructed and film 
level kept steady, even for high current and large amounts of evaporation. 
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2.5 Particle Tracking Velocimetry 

Velocity information is obtained from the e-m cell using a particle tracking method 
(PTV) which is similar to the standard technique known as particle imaging velocime- 
try (PIV). The PIV technique uses a camera to capture images of small particles that 
seed the fluid flow. Two consecutive images are separated in time by a small amount 
At. These images are sectioned into small regions by a discrete grid. Corresponding 
regions from the two consecutive images are compared to obtain the average motion 
Ax of the particles in that region over the time At. Each region is then assigned 
a velocity vector u = Ax/ At, yielding an entire velocity fleld on a grid. There are 
many papers and review articles which describe the PIV technique [11, 22, 23, 24], 
interested readers should refer to these for a complete description. 

Where PIV attempts to track the average displacement of a number of particles 
(usually around 10) in a square region formed by a grid, PTV attempts to determine 
the displacement of individual particles. This allows PTV to obtain flner spatial 
resolution than PIV. Equivalently one could say that PTV has higher vector density. 
Since the algorithms used in PTV and PIV to determine translation are identical 
it might be expected that this increased resolution comes only with increased noise. 
This is not the case, however, because PTV has a built-in self correction that allows 
noise to be suppressed. There is one sacriflce though. Since PTV attempts to track 
individual particles the technique is much more sensitive to particles leaving the 
measurement volume than PIV. This is not a issue in 2D flows such as are studied 
here (except near the boundaries of the images), though in 3D it could be a signiflcant 
problem. A copy of the PTV program is listed in Appendix A. 

There are three main steps to PTV by which one goes from CCD images of 
particles to velocity information: particle identiflcation, neighborhood comparison, 
and matching. Before going into these, the manner in which images are acquired 
should be described. As stated earlier the soap fllm in the e-m cell is seeded with 
small particles. These particles are illuminated by two pulsed Nd:Yag lasers that 
yield 12 mJ of energy per nano-second pulse. Images of the particles illuminated by 
the lasers are obtained using a 30 Hz, 8 bit CCD camera of resolution 768 x 480 
rectangular pixels with aspect ratio 1 : 1.17. The lasers are slaved to the camera so 
that the flrst Nd:Yag laser pulses at the end of the flrst image and the second laser 
pulses a time At later in the second image. In this manner two images of particles 



Chapter 2. The E-M Cell 



22 



Figure 2.8: Timing of the CCD camera frames and laser pulses used in PTV. Frame 1 
and Frame 2 denote a single PTV image pair from which velocity fields are extracted. 

separated by time At are obtained on a camera with a fixed time resolution of .03 
ms. This timing is shown in Fig. 2.8. 

To determine the positions of the particles in the flow the background of each 
image must be subtracted off. There are many routines by which one can do this; in 
this thesis a high pass filter is used since the particles are spatially small compared 
to the image size. Once the background is subtracted the particles in each image are 
found by an exhaustive nearest neighbor searching algorithms. If the pixel is 
found to be non-zero then this serves as the base of a particle group. The four pixels 
at (z + 1, j), (z — 1, j),(z, j + 1), and — 1) are said to neighbor the base pixel, and 
any of these which are non-zero are added to the group. Any neighbor of a pixel in 
the group is then considered, and if it is non-zero and not already part of the group, 
it is added to the group. This process continues until no pixels are being added to 
the group. The final group of non-zero pixels is called a particle. This process is 
performed until all non-zero pixels in each image are accounted for in a particle. In 
what follows the particles in the first image will be indexed by a and those in the 
second image by b. 

This manner of finding a particle does not distinguish between an individual par- 
ticle and particles that are so close together that they form a continuous image on the 
CCD camera. Since the turbulence in the e-m cell is only mildly compressible, par- 
ticles which start initially very close should stay close over the short period of time. 
At, between laser flashes. The indistinguishability, then, should not be an issue. 

There are three quantities which need to be determined for each particle: it's 
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centroid, pixel centroid and size. The centroid is the center of the group of pixels 
which make the particle weighted by the intensity, ijj'^j-), of the pixels in the group (z 
denotes the image). If there are N pixels, indexed by n, in the pixel group of particle 
a in image one, then the centroid, (a;^"\ y^"'^), is given by: 

y-^c ^Uc ) - (1) 

2^n=l\i(n),j{n)) 

Note that where the pixels themselves are discretized on a finite grid, the center of 
mass of a particle need not be if there is more than one pixel contained in it's group. 
This phenomenon is called sub-pixel resolution since it allows particle position to be 
tracked with a resolution finer than the pixel size of the camera. Sub-pixel resolution 
can be used to enhance the dynamic range of PTV measurement, though it is not 
relied upon for results in this thesis. The pixel centroid of particle a is defined as the 
nearest pixel to the particles centroid, and will be denoted as (x^'^\ y^""^) . Finally the 
size of particle a is simply the 2'^'^ moment of the intensity distribution given by 

r(l) \ 



The particle size is used as a filter to discard particles which are too big or too small. 

Once particles have been identified the challenge is to track them from one image 
to the next. This is done by comparing the regions surrounding particles in the first 
image with regions around particles in the second to establish how well they correlate. 
For particle a in the first image and b in the second image the correlation number Ca,b 
is determined by 

?(2) 

^ l^m=-l l^n=-l -'((x(")+m),{j/(")+n)) ({x('')+m),(y('')+n)) 

~(z) 

In the above 21 is called the correlation box size. The intensities, I-j , used in deter- 

(z) 

mining the correlation number are the intensities of the images I^j less their average 
over their respective correlation boxes so that Ca^b assumes a value between —1 and 
1. If the correlation number is close to 1 then the region around particle a is similar 
to the region around particle b. The closer to 1 the more similar the regions. 

One need not compare all particles in the first image to all particles in the second. 
Only a subset of particles within a certain distance of one another need be considered. 
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That is if a particle a is found at the first frame and b is at 

the second, their correlation number need only be calculated if {{x^'^^ — x*^^^)^ + (y^"^ — 
y{b)yy/2 ^ where s is some reasonable threshold displacement based on the RMS 
velocity fluctuations and flash spacing. 

Once all the particle comparisons have been performed, one need only match 
particles in the first frame to those in the second. This is done by an iterative mutual 
maxima technique. Look at correlation number Ca,6 and determine if it is above some 
initial threshold cumit- If it is, then look to see if Ca,6 is the maximum correlation 
value for both particle a and particle b. That is make sure Ca,b > Ca,x for any x ^b 
and Ca,fe > Cyf, for any y ^ a. If so then then Ca^b is the mutual maximum for particles 
a and b and they are considered a match. Thus particle a has moved from position 
the first frame to the second. The fact that we have checked 

that not only is particle a the best fit for 6, but that b is best fit by a is the self- 
correction that PTV has that PIV does not (PIV can only check that a is best for 6), 
and the reason that PTV can achieve higher vector density without much sacrifice in 
velocity resolution. 

Once matched these particles and their correlation numbers are removed from 
consideration. This is done for all correlation numbers above cumiti and all mutual 
maxima are obtained in this way. The cumit is then lowered slightly and the procedure 
performed over all particles which have not already been matched. This is done until 
ciimit hits some specified lowest bound and the particles which have not been matched 
at this point are discarded. 

This leaves a final list of particles which have been tracked from a point in the 
first image to a point in the second. The average velocity of the particle is determined 
by the motion of it's center of mass. This velocity is assigned to the average particle 
position. This yields a field of average velocities which is the final field from the PTV 
technique. These fields are generally interpolated to a finite grid (binned) so that 
derivatives may be taken. This interpolation can be performed by any number of 
weighted averaging techniques. 

2.6 Cell Operation 

This section describes the procedure that was used to employ the features of the e-m 
cell described above. First a marginally thick film is pulled across the square frame. 
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The plastic edges of the frame are dried to remove any fluid bridges that might short 
the electrical current. A small amount of air is removed from the box until the film 
is just level to the eye. This film is then placed above the magnet array and a mixing 
current, generally around 15 mA, is turned on. The feedback loops to inject fluid and 
pressure balance the film are initiated. After a balance is reached the magnets are 
raised (or lowered) to the desired level, and the current is adjusted until the target 
Urms is reached. Particle images are then acquired at a rate on the order of a few Hz 
until a large number (between 500-1000) image pairs are obtained. These pairs are 
then interrogated using the PTV algorithm to obtain velocity information. 



Chapter 3 



Modeling Flows in the E-M Cell 



A frustrating problem that arises when soap films are used as an experimental system 
for studying 2D fluid dynamics is the lack of direct evidence that these films obey the 
2D incompressible Navier-Stokes equation, 

^+Us— = -— + u^ + Fr\ (3.1) 

dt dxs dxi dxl * ' 

In the above equations Ui is the i**^ component of the fluid's velocity field, p is the 
internal pressure field normalized by fluid density, and Ff^^ represents the i'^'^ compo- 
nent of any external force field (per unit mass) acting on the fluid. The constant v 
is the fluid's kinematic viscosity. Einstein summation convention is used, and will be 
used throughout the thesis unless otherwise noted. These equations are the governing 
equations for the time evolution of the velocity field of an incompressible fluid. If the 
soap film does not obey this equation then it is not behaving as an incompressible 
Navier-Stokes fluid and is therefore useless as a system for standard turbulence in- 
vestigations. The purpose of this chapter is to demonstrate that the e-m cell does 
indeed approximate a 2D Navier-Stokes fluid. 



3.1 Introduction 

There are a number of reasons to be skeptical about soap films behaving as an in- 
compressible Navier-Stokes fluids. Most of the problems arise from the presence of 
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Figure 3.1: Thin film interference fringes demonstrate that the thickness of the soap 
film in the e-m cell is not constant but varies from point to point in the flow. 

thickness fluctuations which are caused by the motion of the film. Such thickness 
fluctuations can be imaged using thin film interference as in Fig. 3.1 and are indica- 
tive of compressibility in the soap film. Compressibility constitutes a failure of Eq. 
3.2 in soap films. This has already been discussed somewhat in chapter 2 where it 
was implied that keeping the Mach number small eliminates this problem. This is 
a little misleading; a small Mach number does not mean that there are no thickness 
fluctuations. Rather it means that the thickness fluctuations do not develop shock 
fronts and vary in a smooth manner from point to point in the flow. In the range 
of Mach number present in the e-m cell the thickness fluctuations tend to be around 
20% the mean thickness of the film [11]. One would like to determine if this is small 
enough to allow Eq. 3.2 to hold approximately. 

Aside from the incompressibility issue, thickness fluctuations could cause soap 
films to deviate from a Navier- Stokes fluid in a more sinister way. Recall from the 
discussion in chapter 2 that the kinematic viscosity of the soap film, z/, is dependent 
on thickness. Since the film thickness varies from point to point in the flow, so should 
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the effective viscosity of the fluid. These thickness fluctuations respond to velocity 
gradients in the film, therefore the viscosity is dependent on the local shear rate. A 
fluid with such a shear dependent viscosity is said to be non-Newtonian and does not 
obey Eq. 3.1. Here again one would like to determine if the viscosity fluctuations are 
small enough to be considered negligible and Eq. 3.1 to hold approximately. 

Another problem when dealing with soap films is the external frictional coupling to 
the air. It's presence is important to attain an energy balance in 2D forced turbulence, 
as discussed in chapter 1. Indeed its strength and form should dictate the outer scale 
of the turbulence and affect energy transfer at large scales. Because of its importance 
the effects of this coupling must be modeled and tested. The simplest model of the 
frictional effects of the air is to assume it acts as a linear drag on the film. Therefore 
the external force field F/^* acting on the e-m cell can be broken into two parts, the 
Lorentz force caused by current and magnetic field, Fj, and the air drag F""' = — auj. 
The constant a represents the strength of the frictional coupling of the air to the film. 
This model must be tested if it is to be used in later investigations. 

Though a direct test of the Navier- Stokes equations by inverse methods is not 
easily performed, it is possible to test an equation known as the Karman-Howarth 
relationship. This relationship can be derived from the Navier-Stokes equation with 
a single assumption, and easily tested with data from the e-m cell. It's failure or 
success in describing data from the e-m cell would then constitute an indirect test 
of the Navier-Stokes equation as well as the linear drag model proposed for the air 
friction. 

3.2 The Karman-Howarth Relationship 

Although the derivation of the Karman-Howarth relationship can be found in a num- 
ber of texts on turbulence, it is performed here for two purposes: the relationship is 
used extensively in later chapters and to present notation which will be used through- 
out the thesis. It is also performed with the inclusion of a linear damping term in 
the Navier-Stokes equation to represent the air friction as discussed earlier. The 
derivation here, with the exception of the air drag, closely follows that found in Hinze 
[25]. 

The Karman-Howarth relationship governs the time evolution of the two-point 
velocity correlation, {ui{x)uj{x')) , for a fluid in a state of homogenous turbulence. The 
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brackets (...) represent an ensemble average. This relationship can be derived from 
the incompressible Navier-Stokes equation using only the assumption of homogeneity 
in the following manner. Multiply Eq. 3.1 which is evaluated for the i*"^ component 
of the velocity field at the point x with the j*'^ component of the field at point x': 

Ou d d (P" 

In the above, field quantities at the point x' are denoted by a / and the linear drag 
model has been explicitly inserted into the equation. The fact that the derivative at 
point X commutes with multiplication by fields evaluated at x' has been used to move 

inside spatial derivatives. Incompressibility has also been used in the second term 
on the left-hand-side to move Us into the derivative. 

Add Eq. 3.3 with the corresponding equation evaluated by multiplying Eq. 3.1 
evaluated for the j**^ component of the velocity field at the point x' with the i*'^ 
component of the field at point x. This allows both velocity terms to be brought into 
the time derivative, 

+'^^(^»^) + + ^i-^* + ~ 2aUiUj-. (3.4) 

A coordinate transformation to relative, ri = x\ — Xj, and absolute, = l/2(x^ + 3^1); 
coordinates can now be performed. Grouping the appropriate terms yields: 



d \ d d 



UiUM 



1 92 Q2 

+ -i^^(m.m •) + 2z/— («,«;.) + u'^F, + UiF'^ - 2auiu'^. (3.5) 

Now an ensemble average is performed. Using the assumption of homogeneity 
eliminates the derivative of averages with respect to absolute position, leaving 
only derivatives with respect to relative position, r^. What remains is called the 
Karman-Howarth relationship, 

^(M^n;.) = -A(n,<n;.-UiU,n;) + ^(|)n;.)-^(|)'n,) 

+2z/— (n,n;.) + (n;.F,) + {mF'^ - 2a{uiu'^). (3.6) 

s 
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In the future the n-term two-point velocity correlation functions will be denoted by 



in) 
ij..., 

Eq. 3.6 is given by b\^-, while the first term in the first derivative on the right is given 



b\j = {uiUj...u'i^...) . Using this notation, the correlation on the left hand side of 



'(2) 

by b^^^ 



This relationship can also be used to derive the energy balance equation for ho- 
mogenous turbulence. Energy balance will be used in what follows as a self consis- 
tency check for data that attempts to fit the Karman-Howarth relationship. Taking 
the limit of Eq. 3.6 as r — > (or equivalently as x x')yields 

d 

Q^iuiUj) = Uij - eij + (uiFj) - 2a{uiUj). (3.7) 
In the above the tensors Uij and e^j are defined as 

^^^^H^iry^)-^(p'^^))^ (3-8) 



r— ♦o dr^ 



eij = \im2iy^^{uiu'^). (3.9) 



Taking half the trace of Eq. 3.7 eliminates the pressure term Uij and yields the energy 
balance relationship 

1 d 

2^"rms = -^^rms + (UsFs) - aU^^,, (3.10) 

where the vorticity, u, is defined as the curl of the velocity field (i.e. w = V x n). 
The first term on the right, i^^^rms = amount of energy changed to heat 

by internal viscous dissipation. The second, {ugFs) = ej„j, is the work done by the 
external force. The final term, au^^^ = eair, is the energy lost to the linear drag. 
Equation 3.10 is the statement that the change in energy in the system is simply the 
amount gained from the external forcing less the amount lost to dissipative effects. 



3.3 Testing Karman-Howarth 
3.3.1 Experimental Considerations 

Recall that the objective of the experiments presented in this section is to demon- 
strate that the dynamics of the e-m cell are governed by the Navier-Stokes equation, 
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Eq. 3.1, with the effects of air friction modeled as a Unear drag. This will be demon- 
strated indirectly by showing that the Karman-Howarth relationship, Eq. 3.6, holds 
for homogenous turbulence in the e-m cell. The number of terms which must be 
measured to check Eq. 3.6 can be simplified by using specific characteristics of the 
e-m cell. The first is the elimination of the time derivative. This term can be ignored 
if the turbulence is in a statistically steady state. Since the e-m cell was designed 
specifically to study steady state turbulence it is easy to maintain energy and enstro- 
phy approximately constant. Elimination of the time derivative in this manner is the 
main reason that testing of the Karman-Howarth relationship is significantly easier 
than directly testing the Navier-Stokes equation. 

Another term which can be dropped is the viscous term, if one considers length 
scales, r, greater than the viscous scale rj, ^ (i^^/ej„j)^/^. For typical values of energy 
injection in the e-m cell = 200/im. Since the particle tracking measurements 
focus on the inverse cascade regime, which occurs at length scales greater than a 
millimeter in the e-m cell, most of the measurement resolution lies well above this 
criteria. Since the viscous term is being ignored these experiments can draw no 
conclusions about how thickness changes may be affecting changes in viscosity. Small 
scale investigations, outside the scope of this thesis, would have to be performed to 
draw conclusions about this effect. 

What remains of the Karman-Howarth relationship after using these assumptions 

is 

- bS!,) - |-K) + = K-^') + (««^') - 2«^g- (3.11) 

Normally the assumption of isotropy would also be made to eliminate the pressure- 
velocity correlations on the left-hand-side. Recall from the discussion in chapter 2 
that this assumption cannot be made in the e-m cell due to unidirectional forcing. 
Therefore to check Eq. 3.11 a pressure- velocity correlation must be measured, indi- 
cating that not only is a velocity field needed for the check but a pressure field as 
well. 

To obtain the pressure field, the divergence operator is applied to Eq. 3.1 and Eq. 
3.2 is used. What is left has the form 

V^p = 2A + V-F (3.12) 

where A = — If the Kolmogorov magnet array is used, the divergence 
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of the electromagnetic force on the right may be dropped (see chapter 2). This allows 
the pressure to be approximated from the velocity field using Fourier techniques^. 
For this reason the Kolmogorov array will be used in these experiments. 

With Kolmogorov forcing and the assumptions above, all the terms in Eq. 3.11 
may be measured and tested as an indirect test of the Navier-Stokes equation and 
linear drag model. One final simplification can be made. An exact measure of the 
external electromagnetic forcing is difficult at best. However, since the forcing is uni- 
directional, along the x direction, the force- velocity correlation terms can be dropped 
if the {i-ij) = {y^y) component of the tensor equation is considered. This is easily 
done leaving the final equation to be tested: 

- €ly) - ir(p<) + ir(p\) = -2«C- (3.13) 

All of the terms in the above can be measured, and the constant a can be used as 
a single free fitting parameter. The quality of the fit will determine if the Karman- 
Howarth relationship holds in the e-m cell, and therefore by implication the Navier- 
Stokes equation and linear drag model. Such a detailed comparison between theory 
and experiment has not been performed before for 2D soap film systems. 



3.3.2 The Data 

A single run of the e-m cell using Kolmogorov forcing was performed over a time 
span of ~ 300s during which one thousand vector fields were obtained by PTV. The 
cell was driven at a voltage 40 V with a current of 40 mA oscillating with a square 
wave form at 5 Hz. The magnet array was placed approximately 1 mm below the 
film. This resulted in velocity fluctuations of around 11 cm/s over the time of the 
experiment. A typical velocity field that is obtained by binning the particle tracks is 
shown with the associated pressure field derived using the method described above 
in Fig. 3.2. 

The first order of business is to check that the assumptions of incompressibility 
and that the system is in a steady state are accurate. Shown in Fig. 3.3a is the time 
dependence of the velocity and vorticity fluctuations for the run. The fluctuations are 

^This approximation of the pressure field assumes periodic boundary conditions. Though the 
velocity fields extracted from the e-m cell do not satisfy this boundary condition it is hoped that the 
solution for the pressure field will be insensitive to this approximation away from the boundaries 
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Figure 3.2: Typical velocity (a) and pressure (b) fields obtained from the e-m cell. In 
the pressure field green denotes positive and blue negative values. 
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not exactly constant, but the change is neghgible due to the fact that it happens over 
a long time, i.e. the average time derivative is approximately zero. The steady state 
assumption is therefore approximately correct. The incompressibility assumption can 
be tested by measuring the divergence of the flow, D = V ■ u, and normalizing its 
square by the enstrophy Q = uj'^.^g. This forms a dimensionless quantity which must 
be small if the divergence effect is to be ignorable. In the e-m cell D'^/Q ^0.1 over 
the time of the run as shown in Fig. 3.3b. This indicates that the divergence is not 
overly large and incompressibility can be assumed. Coincidentally this number is also 
close to the Mach number of the system, which the reader will recall was kept small 
for the purpose of minimizing compressibility. 

The final assumption necessary to check before Karman-Howarth is applicable to 
the system is homogeneity. That is the average quantities in the turbulence should 
be invariant with respect to translation. A crude test of homogeneity is to measure 
the mean, {u(x))n, and RMS fluctuation, {{\u{x)\'^) n)K of the velocity as a function 
of position, where N denotes the number of fields the quantity is averaged over. Both 
should be independent of position for homogeneity to hold. Moreover, since the film 
in the e-m cell does not have a net translation, the mean flow everywhere should 
be identically zero. Figure 3.4a shows the the mean flow for the run after having 
averaged over the thousand images (A^ = 1000). One can see that there still exists 
a small mean. Though at first this is discouraging, it is misleading since a finite 
amount of data will almost never converge exactly to zero. Rather the magnitude 
of the fluctuations in the mean shear should decrease as A^~^/^ if one assumes the 
fluctuations away from the mean are behaving as a centered Gaussian variable. To 
this end, the RMS fluctuations of the mean flow, {u{x)) n rms, is plotted as a function 
TV in Fig. 3.4b. It is clear that the magnitude of the fluctuations in the mean is 
decreasing almost perfectly as A^~^/^, indicating that the mean flow as — > oo 
should go to zero as required by homogeneity. 

Although the mean flow is constant (since it's zero) and satisfies the requirement 
for homogeneity, the RMS fluctuations, {{\u{x)\'^) , does not. This can be seen by 
looking at Fig. 3.5 which shows the RMS fluctuations of the two velocity components 
averaged over the thousand images. Although the y component of the velocity fluctu- 
ations is approximately constant, the x fluctuations display strong striations. These 
striations reflect the Kolmogorov forcing, as they must if the electromagnetic force 
is to inject energy into the system. That is, some part of must be non-random 
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Figure 3.3: (a) Time dependence of Urms and Urms for a single run in the e-m cell. This 
demonstrates that the e-m cell is in an approximately steady state, (b) Time depen- 
dence of the enstrophy normalized mean square divergence, D^/u^^g, for a single run 
in the e-m cell. The fact that D'^/u^^g is small indicates negligible compressibility. 
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Figure 3.4: (a) The mean flow in the e-m cell averaged over 1000 vector fields. The 
length of the reference vector in the upper right corresponds to 2 cm/s. (b) The 
decay of the fluctuations in the mean flow as the number of flelds, N, in the average 
increases. The line corresponds to the expected N~^/'^ decay of a centered Gaussian 
variable. 
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Figure 3.5: The RMS fluctuations of (a) and (b) function of position in 

the e-m ceU. Green denotes large values of the fluctuations while blue denotes small 
values. 
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and in phase with the forcing otherwise {F ■ u) = and the e-m cell could not be 
maintained in an energetically steady state. Since the forcing is oscillating in time so 
must this in-phase component; thus it shows up in the RMS fluctuations as a function 
of position and not in the mean flow. Fortunately the magnitude of the in phase part 
of the fluctuations is small, around 1.5 cm/s, compared to the total RMS fluctuations 
of 12 cm/s. Therefore they will be assumed to be ignorable. Other than these os- 
cillations, the cell appears approximately homogenous in the RMS fluctuations as a 
function of position. The assumption of homogeneity can be said to weakly hold for 
the turbulence in the e-m cell. This approximation will be refined in later chapters. 

The Karman-Howarth relationship is now in a position to be tested. For simplicity, 
define 

A. = ^(&S!i (3-14) 

^M--|-K) + ^(A.), (3.15) 

so that the {y,y) component of the Karman-Howarth relationship, Eq. 3.13, may 
be written Ay^y + By^y = —2ab^^y. The three separate terms Ay^y,By^y and b^^y were 
measured and a least squares algorithm used to obtain the a value which best fit 
the measured data to the Karman-Howarth equation. In this case a ~ 0.7 Hz. The 
results of these measurements are shown in Fig. 3.6 a,b,d. In Fig. 3.6c the sum 

Ayy + Byy IS ShOWU. 

First note that By^y, the term involving pressure velocity correlations is non-zero, 
as it would be if the turbulence were anisotropic. This confirms what was earlier 
assumed to be the case, that the unidirectional forcing does not allow isotropy to 
be assumed. Next note that the images in Fig. 3.6c and d have very similar forms, 
namely a central negative trough with two positive peaks on the axis. This is 
evidence that the Karman- Howarth relationship is indeed holding in the e-m cell. 
To get a better feel for the degree to which there is agreement, several cross-sections 
of plots c and d are shown in Fig. 3.7. The noise in the terms Ay y and By y arises 
from the fact that these terms are derivatives, which are always noisy and converge 
slowly in experiment. In spite of the noise there is clearly agreement, and it is therefore 
concluded that Karman-Howarth, and hence the Navier-Stokes equation with a linear 
drag, holds for the e-m cell. 
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Figure 3.6: Measured values of (a) Ay^y, (b) Byy, (c) Ay^y + By^y, and (d) —2ah^y}y 
from Eq. 3.13. Green denotes positive values and blue denotes negative values. 
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Figure 3.7: Cross sections of Ay^y + By^y (— ■) and —2ah^yJ^ (— ) along the lines (a) 
r = r:^ {vy = 0), (b) r = Vy {r^ = 0) and (c) r = = Vy. 
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A quick check to see if this conclusion is correct is to see if the measured coefficient 
for the hnear drag, a, is viable. Using the discussion in chapter 2 the linear drag 
coefficient can be approximated as a = rj/phd. Given a 50 fim thick film, a magnet- 
film distance of 1 mm the coefficient assumes the value 0.36 Hz. This value is at least 
the right order of magnitude, and the difference between this predicted value and the 
measured one may be accounted for by recalling that the drag on the top surface of 
the film has been ignored(see chapter 4). 

3.3.3 Consistency Check: Energy Balance 

The previous section has checked that the (y, y) component of the Karman-Howarth 
equation is consistent with measurements made in the e-m cell. However, the fit to 
the data was somewhat noisy in spite of a thousand fields being used in the average. 
What is needed to bolster confidence in the equations is some sort of consistency 
check. This is provided by the energy balance relationship, Eq. 3.10. 

The energy balance statement for the time independent fiow simply states that the 
energy injected into the system by the electromagnetic force, e^j, must be balanced 
by the energy lost to the air friction, eair, and viscosity, e^. The later two of these can 
now be measured using the definitions of the various e's given earlier, the extracted 
value of a ~ 0.7 Hz and the kinematic viscosity of z/ ~ 0.016 cm^/s. The energy 
dissipated by air is found to be eair ~ 85 cm^/s^, while the energy dissipated by 
viscous forces is ^ 55 cm^/s^. Using energy balance this suggests that the energy 
injected by the electromagnetic force should be ej„j ~ 140 cm^/s^. An independent 
measure of ej„j would then yield a consistency check of the measured value of a and 
the quality of agreement of data to the {y, y) component of the Karman-Howarth 
relationship. 

This check is provided by the (x, x) component of the Karman-Howarth relation- 
ship which allows a measure of emj- This component of the Karman-Howarth equation 
has the form 

A.,. + 5.,. = + - 2a6g) . (3.16) 

Recall that the Kolmogorov forcing is invariant to translation in the ct direction, that 
is, along the forcing. Let us then restrict the displacement vector r to lie along this 
direction. Then (u'^F^) = {u'^F'J) which by homogeneity equals {u^F^) = {u^F'^). But 
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Figure 3.8: Cross section of A^^x + Bx^x {—■) and —2aU^l. i~) along tlie line r = r^. 
(r, = 0). 

^inj = {uxFx). Using these relationships in Eq. 3.16 yields Ax^x + Bx^x = 2ej„j — 2a6^^^ 
along the r = r^, cross-section. Thus the plots of Ax^x + and 2a6^^], should look 
similar to the ones shown earlier except offset by a constant which is equal to 2einj- 
Figure 3.8 shows the plot of the r^. cross-section of the {x, x) components of Ax^x + 
Bx^x and —2ah'^^\. Clearly the plots are offset by positive constant of 2ej„j ~ 240 
cm^/s^. Thus einj ~ 120 cm^/s^ which is close to the expected value of 140 cm^/s'^. 
Systematic data discussed in the next chapter will show that the extracted values of a 
and einj fluctuate around ±20%, thus the measured value of is within experimental 
error of the expected value. This is further supporting evidence that the Karman- 
Howarth relationship does indeed work for data extracted from the e-m cell. 



Chapter 4 

Energy Distribution and Energy 
Flow 

In the last chapter the energy balance relationship, — — eair = 0, was used as 
a consistency check to determine if the Karman-Howarth relationship was applicable 
to data from the e-m cell. Though it was not discussed there, these measurements 
are the first indication that an inverse cascade is present. Recall from the discussion 
in chapter 1 that if an inverse cascade exists then energy is moved from small length 
scales to large ones, away from length scales at which viscous dissipation is effective. 
Thus the bulk of the energy is expected to be dissipated by the external dissipation 
mechanism acting at large scales, which in the case of the e-m cell is the air dissipation. 
This is exactly the result that the measured energy rates in the e-m cell demonstrate, 
i.e. eair > This chapter begins the investigation of the inverse energy cascade in 
the e-m cell by quantifying the length scales over which it exists (it's range), measuring 
how the energy is distributed over these length scales and determining the rate at 
which energy flows through these length scales. 



4.1 Distribution of Energy and the Outer Scale 

The energy spectrum, U{k), provides a means for describing the manner in which 
turbulent kinetic energy is distributed over different wavenumbers (inverse length 
scales) in the e-m cell. It is defined as the average square modulus of the Fourier 
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transform of the velocity fluctuations, U{k) = {u{k)u^{k)y, and it's circular inte- 
gral, E{k) = Jo \k\U{k)de, denotes the average amount of kinetic energy stored in 
wavenumbers of modulus k. If an inverse cascade is present in the e-m cell, the 
expectation is that energy would build up in wavenumbers smaller than the energy 
injection wavenumber kinj. 

This expectation proves to be the case in the e-m cell. Fig. 4.1 are plots of U{k) 
calculated from transforms of the PTV velocity fields for various types of driving 
force in the e-m cell. Note that the peaks corresponding to the injection wavenumber 
differ as expected for the different types of forcing. For example the two square arrays 
produce peaks along the line k^^ = ±ky, with the distance from being dictated by 
the size of the magnets used in the array. Also note that energy contained in small 
wavenumbers is greater than that contained in large wavenumbers, as expected for 
an inverse cascade. This can be better seen in Fig. 4.2 where the circular integrals 
have been taken and a build up of energy at wavenumbers smaller than kinj can be 
seen. 

The Kraichnan prediction for the inverse energy cascade range is that E{k) ~ 
g2/3^-5/3 ^Qj. ^ ^ g^-^^ g typical energy rate (in the case of the e-m cell this is 
eair) [!]• This result is consistent with dimensional predictions for the scaling behavior 
of E[k). Lines corresponding to this prediction have been drawn on the plots in Fig. 
4.2. These lines should not be interpreted as a fit to the data and are drawn only as 
a guide. Only the Kolmogorov data set, (a), appears to be in qualitative agreement 
with the dimensional prediction over slightly less than half a decade of wavenumbers 
below the injection wavenumber. All the remaining data sets have a small range 
directly below the injection wavenumber which could be interpreted as k~^^^, but 
this is quickly lost to a broad peak in the spectrum at small k. This type of behavior 
in the spectrum is in agreement with results reported in [13] where the build up of 
energy at small k is associated with the saturation of energy in the largest length 
scales in the system. Such saturation is not included in the Kraichnan prediction and 
is therefore not indicative of failure of the theory, but rather a failure of the system 
to satisfy the assumptions of the theory^. Since case (a) satisfies the assumptions of 

denotes a complex conjugate. 
■^The assumption which is violated in a saturated system is that the velocity fluctuations are 
homogenous. A saturated system occurs when the outer scale which is determined by external 
dissipation exceeds the system size, as it did for the data in Fig. 4.2(b)-(d). When this happens 
large structures attempt to pack into a small area near the center of the system away from system 
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Figure 4.1: The energy spectrum, U{k), for (a) Kolmogorov forcing, (b) square forcing 
using 6 mm round magnets, (c) square forcing using 3 mm round magnets and (d) 
stretched hexagonal forcing. Green denotes large values of U{k) while blue denotes 
small values. 
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Figure 4.2: The circularly integrated energy spectrum, E{k), for (a) Kolmogorov 
forcing, (b) square forcing using 6 mm round magnets, (c) square forcing using 3 mm 
round magnets and (d) stretched hexagonal forcing. The dashed lines correspond to 
the Kraichnan prediction that E{k) oc k~^^^ [1]. The arrows indicate the injection 
wavenumber 
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Kraichnan's theory one might conclude that the k~^/^ prediction is correct. However, 
it should also be noted that the behavior in the measured E{k) depends on the type 
of window function used when Fourier transforming the velocity fields. In the data 
shown in the Fig. 4.2, a three term Blackman-Harris window has been used. By 
changing this window one could get up to a 20% change in the slope of the energy 
spectra. Due to the limitations imposed by windowing no conclusion may be drawn 
from this data about possible corrections to the Kraichnan scaling prediction. 

The low k limit of the inverse cascade range, denoted kout-, is determined by the 
size of the largest vortices which result from the inverse cascade, and corresponds to 
the low k peak in E{k). The position of this peak should depend on the strength 
of the air friction since this is the large scale external dissipation mechanism in the 
e-m cell. In the chapter 3 the linear damping model for the air friction possessed a 
coefficient a which determines it's strength. Using dimensional analysis, a, and ej„j, a 
length scale called the outer scale can be calculated by Vout = i^inj/d^Y^'^- The outer 
scale represents the size of the vortices at which more energy is lost to air friction 
than is transferred to the next size larger vortices. The outer scale should be related 
to the low k peak in E{k) by kout = '^'n'/fout- 

To check this dimensional prediction, a systematic set of data using the Kol- 
mogorov forcing with various magnet-film distances was taken holding Urms approxi- 
mately constant. Kolmogorov forcing was used so that a and ej„j could be extracted 
using the techniques of chapter 3. Between 400 and 500 vector fields where obtained 
for each magnet-film distance. As in chapter 3, the energy in the e-m cell remained 
approximately constant during the data acquisition time for each run. Table 4.1 lists 
the various constants associated with each of the different data sets. 

The first four data sets listed in Table 4.1 may be compared for the purpose of 
error analysis. The first and second of these were obtained using identical values of 
the external control parameters and thus the extracted values of a and emj should 
be identical. This is found to be the case up to two significant figures for the first 
two data sets. The third and fourth data sets are also taken under identical control 
conditions different from those used for the first and second data sets (the magnet-film 
distance was slightly smaller). In these data sets the values of a and emj vary around 
the mean by about 10%. Thus, to be conservative, the error in the two quantities a 

boundaries. This will be investigated in some detail later in this chapter. 
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etnj (cmVs^) 


a (s-^) 


Urms (cmVs^) 




2-71 / kout{cm) 


rint{cm) 


case 


63 


0.45 


7.81 


2049 


3.14 


0.63 


a 


63 


0.45 


7.80 


2014 


3.14 


0.64 




101 


0.6 


8.96 


2901 


2.72 


0.59 


b 


110 


0.65 


9.18 


3103 


2.63 


0.58 




150 


0.9 


9.18 


3579 


1.56 


0.51 


c 


154 


1.25 


7.81 


3211 


1.31 


0.41 




197 


1.55 


7.95 


3567 


1.31 


0.38 


d 



Table 4.1: Global constants for several runs of the e-m cell using Kolmogorov forcing 



and einj will be assumed to be as much as 20% of the measured value. 

The E{k) measured from the data sets labeled a,b,c and d are displayed in Fig 
4.3. The low wavenumber peak, kouti moves to smaller and smaller wavenumber as 
a decreases. This is in qualitative agreement with the dimensional prediction. Using 
the position of this peak the outer scale is calculated by Tout = ^-k jk^t and compared 
to that obtained from the dimensional prediction using the measured values of a 
and einj- The results of this comparison are shown for all the data sets in Fig. 4.4. 
The vertical error bars reflect the propagated error of e^j and a while the horizontal 
represent the discretization inherent in finite Fourier transforms. From Fig. 4.4 one 
can see that the dimensional prediction of the outer scale is not inconsistent with the 
measured outer scale using the low k peak in E{k\ 

The measurements presented here clearly indicate that an inverse cascade is 
present in the e-m cell for all types of forcing. The inverse cascade range is shown 
to exist over wavenumbers k such that /cout < k < k^j, where kinj is the energy in- 
jection wavenumber determined by the electromagnetic forcing and kout is the outer 
wavenumber determined by the external dissipation, kout is found to be not inconsis- 
tent with the dimensional prediction, kout ~ i^inj/ct^Y^'^- No strong conclusion can 
be draw about the manner in which energy is distributed over this range due to win- 
dowing difficulties, however the Kraichnan prediction of E{k) ~ k~^^^ superficially 
holds for data sets in agreement with the assumptions of the prediction. 

Before leaving this section, the systematic data set used in obtaining the outer 
scale allow the external dissipation to be compared to the predicted value a = 
Vair/phd. The measured values of a versus the magnet-film distance are shown in 
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Figure 4.3: The circularly integrated energy spectrum, E{k), for the four cases of 
Kolmogorov flow labeled in Table 4.1. 



Figure 4.4: Comparison of the outer scale obtained from the energy spectra by rout = 
'^'^/kout with that obtained using the dimensional prediction rout = i^inj/ct'^Y^'^ for all 
of the data sets in Table 4.1. 
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Figure 4.5: The measured linear drag coefficient, a, versus the magnet-film distance, 
(i, for the data sets reported in Table 4.1. The dotted line represents the fit a = 
Vair/phd + C with C = 0.25 Hz. 

Fig. 4.5. The vertical error bars again refiect the 20% error in measurement while 
the horizontal error bars denote the limit of control over magnet-film distance. A line 
corresponding to rjair/ phd + C, where C = 0.25 Hz is also plotted with the data and 
for the most part is within error of the measured values. According to this data the 
prediction for the magnitude of the linear dissipation must be offset by a small posi- 
tive constant to be accurate. This constant can be accounted for by recalling that the 
effect of air friction on the top surface of the film has been ignored. Approximations 
of the frictional force on the top surface of the film indicate that the measured value 
of the offset is appropriate, though more experimentation needs to be done to more 
accurately account for this offset. 



4.2 Energy Flow 

A simplified viewpoint of how kinetic energy might be transfered through the inverse 
cascade range in the e-m cell is to imagine that the energy is fiowing like a liquid 
through a pipe. Energy produced by the forcing is poured in at one end of the pipe 
which characterizes the injection scale. It is then moved along the pipe to larger 
length scales by the mixing of fiuid (i.e. by vortex cannibalization) . Finally energy 
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is exhausted from the pipe at the opposite end which characterizes the outer length 
scale. Since the total energy in the system is constant, the amount of energy poured 
into the pipe must be equal to the amount exhausted from the pipe. It is also expected 
that the rate of energy being poured into the pipe is equivalent to the rate of energy 
transferred across any length scale in the middle of the pipe. That is to say that the 
energy flux through the pipe is not dependent on the position in the pipe and remains 
constant. 

This viewpoint is a severe simplification of what actually happens in 2D turbu- 
lence. First it ignores losses of energy to viscous dissipation and assumes all energy 
is lost to external dissipation. Since viscous losses are presumably smaller than losses 
to external dissipation, let us accept for now that this simplification is valid. A more 
important simplification, the one which is of concern in this section, is that the pipe 
doesn't leak. That is there are no holes drilled along the length of the pipe. That 
is to say energy cannot be exhausted from the pipe by external dissipation until the 
outer length scale is reached. To hope that this is actually the case in the e-m cell, 
or for that matter any other laboratory 2D turbulence system is quite a stretch. 

It is more likely that there exists a range in the pipe, probably close to the injec- 
tion scale since external dissipation increases at large length scales, over which the 
amount of energy lost to leaks is negligible compared to the energy flux through the 
pipe. This range will be called an "inertial" range since the energy flow is almost 
entirely dictated by the fluid's inertia and not the energy dissipation. The inertial 
range is of interest due to it's universality. Presumably two 2D turbulent systems 
with completely different external dissipation mechanisms will behave identically in 
their inertial ranges. To use the pipe analogy, the inertial range of both systems 
is completely closed and thus fluid mixing and energy flux should be the same in 
this region. Outside of this range the pipe leaks, and the energy flux depends on 
how many holes of what size are drilled at what position in the pipe. Therefore the 
characteristics of the turbulence may not be universal outside of the inertial range. 

To determine if a range is inertial or not, a measurement of energy flow must be 
made. One way in which energy flow may be characterized is by the third moment 
of velocity difference. The third moment, labeled S^^\r) for now, can be thought of 
as the average energy per unit mass advected over a circle of radius r centered at x 
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per unit circumference of the circle: 



S(^\r) = ^ r {E{x + r)u{x + r)-r)de. 



2-Kr Jo 

Assume also that the reference frame is moving so that the velocity u{x) = 0. As- 
suming homogeneity the circle can be placed anywhere in the turbulence, so that 
S^^^ does not depend on x. E in the above is the energy per unit mass. If the third 
moment is positive, then on average energy is being advected from inside the circle to 
outside the circle, and vice versa if the third moment is negative. It is interesting to 
note what happens if r is assumed to be in an inertial range. If this is the case then 
no energy is lost to external dissipation at that length scale. Since the energy held in 
the turbulence at a scale r is in a steady state then all of the energy injected by the 
forcing into the circle must be advected over the surface of the circle. Call e the rate 
of energy injection per unit mass into the system. The rate at which energy is injected 
into the circle is then vrr^e. Replacing the integral with this yields 5'^^-*(r) «i er. Thus 
a linear range in the third moment indicates inertial behavior of the energy transfer. 

The above derivation of a linear behavior in the third moment for an inertial range 
can be put on much more solid foundation. Indeed the third moment is one of the 
few quantities for which an exact prediction can be derived from the Navier-Stokes 
equation. This derivation was first done by Kolmogorov for 3D homogenous and 
isotropic turbulence. To apply to results from the e-m cell Kolmogorov's derivation 
must be relaxed to the case of 2D homogenous but anisotropic turbulence with an 
external linear drag. This relaxation is given below. Following this are measurement 
and analysis of the third moment in the e-m cell for the data sets in Table 4.1. 

4.2.1 The Anisotropic Third Moment 

The starting point for the derivation of the third moment relationship for homogenous 
anisotropic 2D turbulence is the Karman-Howarth relationship, Eq. 3.6, which was 
derived in section 3.2. All of the notation and conventions used in that section will 
be carried over without alteration. The notation and the first step of the derivation 
follows that given in a recent paper by Lindborg [26] . 

Add Eq. 3.6 evaluated for (i, j) to that evaluated for (j, i), and use Eq. 3.7 to 
obtain 

2n--2e^'^) = —B^^H-B^^-—P-—P- 
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7^^(«-0 = 7^^S- (4-2) 



In the above, the n term moments Blj\{r) = {6ui6uj...Suk) of velocity difference 
6ui = u[ — Ui have been defined. Also defined are Pi = {uip') — {u'-p) and Wij = 
{6ui6Fj) with 6Fj = Fj — Fj. Note that the homogeneity assumption has been used 
in a number of places in the above step. Most notably, it sets 

drs^'''' ''''' drs 
Contracting Eq. 4.1 with the unit vectors nj(= r^/r = Ti/|r|) and Uj and using 

the identities 

n,n,Blf = (4.3) 

= ^(^^^.^S) - IbHII (4.4) 
d d 1 



(4.6) 

where the subscripts r and t denote longitudinal, i.e. along r, and transverse direc- 
tional coordinates, we obtain 

-u(-( _ R(2)^ ,1^ o(2) o(2)\ 



+«i?(^)_ln,n,(PF., + PF,,). (4.7) 

Eq. 4.7 is now in a form which may be easily integrated over a circle of radius 
r. This procedure, along with incompressibility and the assumption of homogeneity 
eliminates the pressure terms. Using the divergence theorem and rearranging the 
terms yields 



r Jo r Jo ot r Jo 

+2u (-Sif + ^Sif + - rdr'{Sil\r') - S^\r'))/r' 
\r or r Jo 

— I 

2iTr Jo 



f — r dr' r'ddninj{Wij{r') + Wji{f')). (4.8) 
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Here the circular averages of velocity moments have been denoted as S^jf^^r) = 
^ Jq^ rd6Bl^\{r). Though this equation seems to explicitly contain an er term, 
it is somewhat superficial as this term exactly cancels with terms contained in the 
first and final expressions on the right hand side. Removing these terms yields 

Sl'lir)-'^ J^dr'sShr') = dr' f\'den,n,{u[F, + u,F'. + u'^F, + u,Fl) 

+2v + ^^Sl^ + \ [dr'iS^Xr') - (r'))/r') 

+ f-l: + — 1 r^^' rr'dehfli?). (4.9) 
\Txr ot vrr J Jo Jo 

Note that using the notation introduced here b^^} (r) is simply the two-point longitudi- 
nal velocity correlation. The terms on the left hand side are defined as the anisotropic 
third moment of velocity difference, Sj^\ and up to a constant play the role of the 
third moment of the longitudinal velocity difference in the fully developed isotropic- 
homogeneous turbulence derived by Kolmogorov [9] . The terms on the right account 
for the energy flux at some length scale due to external forces. Eq. 4.9 is essentially 
the scale-by- scale energy balance relationship for the system. For large r the viscous 
term in Eq. 4.9 can be ignored, as can the time derivative if the system is in an 
energetically steady state, leaving the flnal form which will be used in this thesis 

Slllir) - - r dr'sSl (r') = dr' ^ r'dernnj {u[Fj + u,F' + u\F, + UjF[) 

r Jo 2vrr Jo Jo j j 

+ — / dr' / r'deh^^lir'). (4.10) 
7rr JO Jo 

Limits of this equation are now ready to be taken to establish that a linear range 
can exist. First note that the force- velocity correlation term in Eq. 4.10 (flrst term on 
the right) should decay in magnitude as 1/r for r ^ rj„j for F periodic in space. One 
limit which can be considered is the case where this force-velocity term is neglected. 
In this case the only remaining term arises from the linear dissipation. Assuming 
S^{r) oc r^/^ over some range of length scales for r > rj„j, then h^^}{r) = u^^g — 
S^{r)/2 oc ul^g — Ar"^/^ in that range, where A is a constant. Since the longitudinal 
velocity correlation, h^} (r), remains positive the flrst term is dominant. Thus h^^} (r) ^ 
ulrns ^his range. Inserting this approximation into Eq. 4.10 and integrating yields 
S^^ = 2au1^gr = 2eairr, which is the extension of the earlier mentioned Kolmogorov 
4/5 result for 3D turbulence to 2D anisotropic turbulence. Another limit of interest is 
when both the force- velocity and velocity-velocity correlations have disappeared. In 
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this case the integrals on the right hand side become constant and the third moment 
decays as 5*^^^ ^ r~^. 

Notice that the e found in the third moment relationship is eair, the energy dissi- 
pated by the e-m cell's external dissipation mechanism of air friction. One might have 
expected that this should be e^j, the total energy injection rate. Recall that in the 
previous discussion of the third moment, energy lost to viscous forces were ignored so 
that einj = eair- When viscosity is reintroduced, then the energy flowing over a circle 
of radius r is the energy injected less the amount dissipated by viscosity, i.e. e^j — e^. 
By energy conservation this is just eair, which dictates the remaining energy which 
must flow over the circle. This is why eair determines the third moment and not einj. 

4.2.2 Homogeneity 

Before sliding headlong into Sj^^ measurements and blindly searching for positive 
linear ranges, a word of caution is warranted. The assumption of homogeneity has 
been used a number of times in the preceding derivation of Sj;j^\ not to mention the 
fact that it was used in deriving the Karman-Howarth equation. This makes S^^ 
measurements extremely sensitive to inhomogeneity in the e-m cell.^ 

Fortunately the derivation of S^^ has provided a simple test of homogeneity. If 
the turbulence in the e-m cell is homogenous then the following should be equivalent 
representations for the anisotropic third moment, S^^: 



The latter two forms will be denoted J(r) and K{r), respectively. Plots of all three 
of these quantities for the data sets labeled (a)-(d) in Table 4.1 are shown in Fig. 
4.6. It is clear that only case (d) of extremely heavy damping (large a) produces 
approximate agreement for all three forms at length scales larger than the injection 

■'There might be some concern that results presented in chapter 3 might be inaccurate due to 
inhomogeneity. RccaU that in that chapter homogeneity was assumed and not exactly checked. 
Using the analysis presented in this section on the data of chapter 3 demonstrates that these data 
sets are approximately homogenous. 
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scale. It is therefore the only approximately homogenous data set. It is also evident 
that case (a) and (b) are strongly inhomogeneous, and case (c) is marginal. 

To understand the nature of the discrepancies, ensemble averages of individual 
velocity fields were performed over images. Again the ensemble averaged velocity 
at any given point in the flow, {u{x))n, though not identically zero was found to 
decrease in magnitude as N~^^'^ for all of the data sets indicating negligible mean 
flow in the system. The inhomogeneity, then, stems from the spatial variation of the 
velocity fluctuations, {{\u{x)\'^) pf)^ , which is shown in Fig. 4.7 for the four data sets 
(a)-(d). The oscillating light and dark bands, corresponding to high and low u^ms 
and inhomogeneity at the injection scale are again present in all four sets to a greater 
or lesser extent. As before these oscillations will be ignored. A closer inspection of 
the Urms fields for the four data sets also reveals a large-scale inhomogeneity which 
increases in magnitude as a decreases. Note that for the most weakly damped case 
(a), the fluctuations near the corners are weak compared to those near the box center. 
This large-scale inhomogeneity is the main source of the discrepancies for the three 
different forms of 

The reason that the velocity fluctuations begin to form this large scale inhomo- 
geneity for weak damping is connected to the growth of the outer scale of the turbu- 
lence. As discussed before the outer scale indicates the largest sized vortices present 
in the turbulence. This can be seen in Fig. 4.8 which shows typical streamlines for the 
four cases (a)-(d). Clearly the heavy damping, case (d), has many small vortices but 
few large ones compared to the case of weak damping, (a). These large vortices prefer 
to exist in regions removed from solid boundaries, otherwise a large shear builds up 
between the vortex and the boundary and quickly dissipates the vortex. Since the 
largest vortices are of diameter Vout, this preference causes an absence of fluctuations 
for distances smaller than rout from the wall. Should this boundary region invade the 
PTV measurement homogeneity is sacrificed. 

For all of the data sets in Table 4.1 the distance between the PTV measurement 
area and the boundary was approximately 1.5 cm, this sets the largest outer scale 
possible before sacrificing homogeneity. Table 4.1 reveals that the spectrally measured 
outer scale, 27r /kout, of the case (a) and (b) are well above this value explaining their 
strong inhomogeneity. Case (d) is below this value therefore the measurement volume 
is homogenous. Case (c) has an outer scale just exceeding this limit, which explains 
it's marginal homogeneity. 
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Figure 4.6: Comparison of S^\r) (o),J(r) (<), and K(r) (o) for the data sets labeled 
(a)-(d) in Table 4.1. 
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Figure 4.7: Spatial variation of velocity fluctuation for the four data sets labeled in 
(a)-(d) Table 4.1. Green denotes large values of the fluctuations while blue denotes 
small values. 
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Figure 4.8: Typical streamlines for the four cases labeled (a)-(d) in Table 4.1. 
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4.2.3 The Inertial Range and The Integral Scale 

In the last section only one of the data sets which were analyzed strictly satisfied the 
homogeneity condition, case (d). Case (c) was marginal in it's homogeneity, so it can 
be considered as well. Considering only these two cases it is apparent that there is no 
inertial range in the e-m cell. Recall that from Eq. 4.10 an inertial range is indicated 
by a linear range in S'p)(r). Neither case (c) or (d) shows such a range in Sj^^\ Since 
both (c) and (d) have an inverse cascade range, this is the first indication that the 
inertial range is behaving distinctly from the inverse cascade range. However, since 
neither case (c) or (d) has an extensive inverse cascade range, it would be helpful to 
determine if either case (a) or (b), both of which exhibit large inverse cascade ranges, 
has an inertial range in spite of their lack of homogeneity. 

To that end consider Eq. 4.10. The left-hand side of this equation can be rep- 
resented by any of the three forms S^^\r), J{r), or K{r) from the last section if 
the turbulence is homogenous. For the inhomogeneous case it would be useful to 
find if any of these three forms satisfy Eq. 4.10, effectively relaxing the condition 
of homogeneity on the left of the equation. To test this idea the right-hand side of 
Eq. 4.10, denoted R{r), was independently measured and is displayed in Fig. 4.9 
along with K{r) the third representation of the anisotropic third moment for the four 
cases discussed earlier. The strongly and moderately damped cases produce mod- 
erate agreement between K{r) and R{r) , while for the weakly damped cases there 
is some disagreement which increases for large scales. The agreement for all cases, 
however, is better for K{r) than if S^^\r) had been used to represent the left hand 
side. One can weakly conclude then that the homogeneity assumption used to obtain 
Eq. 4.10 can be relaxed on the left for cases of moderate inhomogeneity if the K{r) 
representation of the anisotropic third moment is used. For this reason we call K{r) 
the quasi-homogenous part of the third moment. 

Now the determination of the inertial range for weakly inhomogeneous turbulence 
boils down to whether or not a linear range exists in the quasi-homogenous part 
of the anisotropic third moment for r > rj„j. Clearly, even for the cases of weak 
damping, this does not happen. Indeed over the majority of the range displaying 
inverse cascade, the quasi-homogenous part of the anisotropic third moment decays 
as r~^, which is indicative of a linear dissipation dominated regime. Therefore, none of 
the inverse cascade ranges produced in the e-m cell are inertial, they are all dissipation 
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Figure 4.9: Comparison of K{r) (o) with the independently measured right hand side 
of Eq. 4.10 ,R{r) (— ),for the data sets labeled (a)-(d) in Table 4.1. 
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dominated. 

In the previous section, the outer scale of turbulence was determined to behave as 
i^out i^^y^"^- This length scale, obviously, does not determine the upward extent 
of the range over which energy transfer is inertial, otherwise, a few of the data sets 
should display a linear range. It is reasonable to ask what condition is required 
in order for any linearly damped 2D system to have an inertial range. Equation 
4.10 indicates that if a linear range is to exist it arises from the dominance of the 
second integral on the right-hand side over the first integral at large scales. The 
first integral, the force- velocity correlation, decays (as 1/r) at length scales larger 
than the injection scale. Thus, if a typical length scale of the second integral greatly 
exceeds the injection scale, the third moment should have some linear range. The 
second integral contains the longitudinal two-point velocity correlation as one of its 
functional arguments. This suggests that a measure of the typical length scale of this 
integral is the so-called integral scale. 



The size of the integral scale for the data sets is displayed in Table 4.1. Note that all 
but the most weakly damped data sets have integral scales smaller than the injection 
scale [vinj ~ 0.6 cm). For the case of very weak damping the integral scale has just 
exceeded the injection scale. This is reflected in the quasi-homogenous part of the 
third moment for the weakly damped cases (a) and (b), which show an extended 
plateau right after the injection peak. To better visualize this plateau growth, the 
right hand side of Eq. 4.10 for the four cases is displayed in the Fig. 4.10 with the 
plots normalized by the value of the peak after the injection scale. Such a feature is 
absent in the more heavily damped cases. Thus, the integral scale seems to govern 
the upward extent of the inertial range, in contrast to the outer scale that governs 
the upward extent of the inverse cascade in 2D. 

One might believe that the integral scale and the outer scale should be linearly 
proportional to each other. This may not be the case. Figure 4.11 is a comparison of 
the the outer scale calculated using the dimensional argument with the integral scale. 
Though the error in the plot does allow for the possibility of a linear fit, it is more 
likely a power law growth as shown by the rf^^ curve drawn in the figure. Within 
the error indicated on the plots, the integral scale seems to behave as the geometric 
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Figure 4.10: The right hand side of Eq. 4.10 ,-R(r),for the data sets labeled (a)-(d) in 
Table 4.1: (a) o, (b) >, (c) <, (d) o. i?(r) has been normalized so that the peak value 
just after rinj is unity. 
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Figure 4.11: The dimensionally predicted outer scale, (emj/^'^)^^^ vs. integral scale 
Tint (inset is the same plot on log- log scales). The line corresponds to the power law 
fit of rl^. 

average of the outer scale and injection scale (see the inset), 



Tint = y/rinjTout = 3 • (4.14) 



It seems that this relationship should be predictable by inserting finite ranges in the 
energy spectrum and inverse transforming to get the two-point correlation. At this 
time such a calculation has not been performed. 



Chapter 5 

High Order Moments 



Recall from the introduction that the inverse energy cascade range was expected to 
have two important properties: locality and scale invariance. The experiments in the 
e-m cell do not have enough inverse cascade range for any conclusion to be directly 
drawn about either of these properties. However, some indirect conclusions pertaining 
to scale invariance might be drawn if one assumes that certain measured properties 
of the moments of velocity difference can be extended to the case of a large inverse 
cascade range. Displaying these properties and demonstrating how such conclusions 
may be drawn is the purpose of this chapter. First, to facilitate the extension of mea- 
sured results, a brief discussion which describes the expected behavior for moments of 
velocity difference assuming scale invariance is presented. This discussion is similar 
to one in [9]. 

5.1 Scale Invariance and Moments 

In the introduction the velocity difference over a length scale r, 5u{r) = {u{x + "f) — 
u{x)), was said to be scale invariant in the inverse cascade range (Note that though 
Su depends on the absolute position, x, any statistical quantities formed from 6u 
does not if homogeneity is assumed). Here, scale invariance means that there exists 
a unique scaling exponent h such that P{6u{Xf)) = P{X'^6u{r)), where P denotes 
a probability distribution function (PDF). For the moment, consider homogenous 
isotropic 2D turbulence. Instead of using both components of Sul-f) in P, isotropy 
allows for the simplification to only a single component. For various reasons this 
component is usually the longitudinal component, Su\\{r) = 5u{r) ■ f. Also note that 
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isotropy allows dependence on fto be reduced to only dependence on r = |?^. From 
this the scaling relationship becomes P{6u\\{Xr)) = P{X'^u\\{r)). 

The n}^ order moment of longitudinal velocity difference is defined as 

Si±{r) ^ ((5n||(r))") = 1^^ d5u\\{r) {5u\\{r)r P {5u\\{t)) . (5.1) 

J — oo 

Note that these functions were already defined in chapter 4. Since longitudinal 
fluctuations will only be considered here the notation can be simplified by defin- 
ing Sn{r) = S*^"^ (r). The scale invariance of the PDF's directly results in the scale 
invariance of the moments of velocity difference, with the scaling exponent nh for the 
^th QYf^QY moment. That is 

/ + 00 
d5u\i{\r) (5n||(Ar))"P(5n||(Ar)) 
-oo 

/H-oo 
dX^6u\\{r) (A'^5M||(r))"P(A^'5M||(r)) 
-oo 

/ + 00 
X^dduiiir) (5M||(r))"A-^P(5M||(r)) 
-oo 

= A"^5„(r). (5.2) 

Thus scale invariance implies that the moments of longitudinal velocity difference 
behave as Sn{r) oc r"^. 

Recall from the discussion in chapter 4 that there exists an exact result for the 
third moment of velocity difference in an inertial range, S^i^r) = |er. If an inertial 
range exists, and if it is scale invariant, then this exact result fixes the scaling exponent 
h = 1/3 and the final expectation for the scaling behavior of the n}^ order moment 
is that Sfi (r) oc r"". This result was first reached by Kolmogorov in his 1941 theory 
of homogenous isotropic turbulence and will be called the K41 theory. The result 
can be fleshed out a little more if one assumes that the only variables of importance 
in the inertial range are the constant energy flow rate, e, and the length scale, r. 
Under these assumptions dimensional analysis predicts Sn{r) oc C.„(er)"'/^, where the 
dimensionless C„ are assumed to be universal constants. Incompressibility sets Ci = 
and the previously mentioned exact result sets C3 = 3/2. 

This analysis yields a simple way in which the scale invariance of the turbulent 
flelds can be verifled. Simply make sure that the PDF's of longitudinal velocity 
difference for various r in the inertial range, when properly normalized, collapse to 
a single curve. Equivalently, make sure that the n}'^ order moments of longitudinal 
velocity difference scale as r"/^ for all n. It should be pointed out that this analysis 
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holds for homogenous isotropic 3D turbulence, except for a change in coefficients 
Cn- Although one might expect the velocity differences in 3D turbulence to be scale 
invariant as expressed above, it is an experimental fact that the moments steadily 
deviate from the expected scaling behavior of the K41 theory as the order of the 
moment is increased. The possibility of such deviations happening in 2D turbulence 
is still an open experimental question, though what follows in this chapter might be 
considered clues as to what the answer may be. 

5.2 Disclaimer 

The analysis in this chapter must begin by pointing out a couple of reasons not to 
extend certain conclusion presented in it to cases of 2D turbulence beyond the e-m 
cell. The first of these reasons stems from the type of statistical quantities which 
will be used in the determination of scale invariance. The analysis in the previous 
section was phrased solely in terms of longitudinal velocity differences. In order to go 
from full velocity differences to longitudinal differences without loss of information 
it is necessary that the turbulence be fully isotropic. The fact that the forcing is 
unidirectional, as discussed in chapter 2 does not allow isotropy to be assumed in the 
e-m cell. This means that the analysis given above must be done for fully anisotropic 
turbulence to be absolutely correct. 

Unfortunately, this analysis becomes prohibitive for anisotropic turbulence as the 
order of the moments increase since the number of moments containing transverse 
velocity differences that need to be measured grows. Earlier experiments show that 
the deviations from scaling prediction in 3D become apparent only for moments of 
large order. Assuming that 2D might be similar means that a large number of quan- 
tities would have to be accurately measured and compared, making fully anisotropic 
analysis exceedingly difficult. The fact that anisotropy does not allow the dependence 
on the difference vector r to be simplified to dependence on r = r exacerbates this 
difficulty. From an experimental point of view this is devastating since dependence 
only on r allows statistics to be averaged over circles, increasing the number of points 
used in evaluating the PDF's by 27rr for each scale r. Without this buffering of the 
statistics one cannot hope to obtain enough data to measure high order moments of 
the PDF. For these reasons a fully anisotropic analysis must be abandoned. All results 
in this chapter use only longitudinal differences, and it is hoped that the anisotropy 
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will not seriously affect the final results. 

The second reason was hinted at in chapter 4. The results that will be discussed 
will be coming from an inverse cascade range that is not inertial. This means that 
the form of the external damping may be affecting the results. Since the external 
damping in the e-m cell is known to have a linear form, the results can strictly 
be said to only apply to other linearly damped 2D turbulent fluids. This fact, along 
with the e-m cells anisotropy, raises serious questions about the universality of certain 
results which are obtained. In spite of this, the author feels that many of the results 
presented are robust due to similarities with data obtained from other experiments 
[5] and simulation [27]. 

5.3 The PDF of Longitudinal Velocity Difference 

Of the data sets presented in Table 4.1 the only two to display homogeneity or 
marginal homogeneity by the analysis of chapter 4 are cases (c) and (d). These 
will be the main data sets from which conclusions in this chapter will be drawn. In 
particular, focus will be given to (c) since it has the largest inverse cascade range of the 
two as shown in Fig. 4.3. The PDF of longitudinal velocity difference, P = P{6u\\{r)) 
was calculated in a straightforward manner from data set (c) and is presented in the 
color plot of Fig. 5.1. Several cross-sections of the plot for various r are shown in 
Fig. 5.2. 

It is clear from the cross-sections that P is approximately Gaussian at all r. Of 
course it cannot be perfectly Gaussian since there must be some third moment as was 
measured in chapter 4. The magnitude of the odd moments, such as the third, must 
be small compared to the even order moments for P to have such a strongly Gaussian 
character. The tails of P do not seem to strongly deviate from Gaussian decay into 
either exponential or algebraic decay at any r in the inverse cascade range. In 3D, 
deviations of the high order moments from the K41 theory is associated with slower 
than Gaussian decay in the PDF tails. This behavior is commonly termed "intermit- 
tency" since it indicates intermittent bursty behavior in the velocity field. That no 
such deviation in the PDF tails is seen in these experiments is an indication that the 
scale invariant result may hold. The approximately Gaussian PDF's measured here 
are in agreement with both recent experiments [5] and simulations [27]. 

The moments of the longitudinal velocity difference are calculated from P using 
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Figure 5.1: P{6ui{r),r) calculated from data set (c) in Table 4.1. Divisions in the 
coloration increase on an exponential scale. The injection and outer scale are marked 
by lines, in between which is the inverse energy cascade range. 



Figure 5.2: Cross sections at various r for P{6ui{r),r) shown in Fig 5.1. 
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Eq. 5.1. Since the two lowest order moments in 3D do not measurably deviate from 
the K41 theory, and it is expected that this will be the case in 2D, the low order 
moments are analyzed in the e-m cell first. Displayed in Fig: 5.3 are S2{r) and S^i^r) 
for the data set (c). Clearly there is no range in the second moment which scales 
as r^/'^ in between the injection and outer scale as K41 predicts. There is no linear 
behavior of 5*3 (r) in this range either. This later result is hardly surprising considering 
that a linear range is indicative of an inertial range which has already been shown 
not to exist in the e-m cell (see chapter 4). 

These graphs clearly violate the scaling prediction in the K41 theory, therefore 
the e-m cell's inverse cascade range is not scale invariant. It might become scale 
invariant if an inertial range was allowed to build in the cell. To better understand 
what might happen, the higher order moments (i.e. n > 3) of the longitudinal velocity 
difference are evaluated. These moments are more easily displayed if they are made 
dimensionless by dividing out the appropriate power of the second moments. Define 
Tn{r) = S'„(r)/(5'2(r))"/^, so that T3 is the skewness, T4 is the flatness etc. These 
normalized moments are shown in Fig. 5.4 for 4 < n < 11. The error bars are 
calculated by truncating the PDF wherever noise dominates the PDF measurement. 

First consider the moments of even order. With the exception of a blow up at 
small r, the even T„(r) are essentially constant for r > rj„j with the exception of 
logarithmically small fluctuations between rinj and Vout- The constant values that 
the even T„ assume at large r are only slightly less than those of a pure Gaussian 
distribution, namely 

n\ 

" 2W2)(n/2)!' ^^-^^ 

for n > 2. These are shown on the even plots of Fig. 5.4 as dotted lines. This confirms 
the earlier visual observation that P was approximately Gaussian. The large blow 
up at small scales r < r^nj is due to poor statistics. The fluctuations of the even T„ 
at scales in between rj„j and Tout seems to force the value slightly higher than the 
Gaussian value. It should be pointed out that a similar rise is seen in [5] though 
without the fluctuating character. Fluctuations are also seen in [27] near the outer 
scale, though these settle back down to the Gaussian values once the inertial range is 
reached. Since the fluctuations are small the Gaussian values will be assumed to hold 
for all r > Vinj. From the measured constant values of the even T„ the higher order 
even moments in the e-m cell can be approximated as Sn{r) ~ an('S'2(r))"/^ with the 
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Figure 5.3: (a) 52 (r) (log-log) and (b) 5*3 (r) (lin-lin) calculated from data set (c) in 
Table 4.1. 
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Figure 5.4: The normalized high order moments, T„(r), evaluated from data set (c) 
of Table 4.1 for 4 < n < 11. 
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Figure 5.5: Tn{r)/hn for odd n > 3 evaluated using the data set (c) in Table 4.1. 

Qn assuming the values of a Gaussian distribution. 

Now consider the odd moments. Unlike the even moments, the odd T„ display a 
complicated behavior for r > rj„j. However, up to a multiplicative constant, which 
will be denoted 6„, the odd T„ have similar behaviors. To see this the multiplicative 
constant has been removed and all odd T„ for n > 3 have been plotted in Fig. 5.5. The 
hn where chosen to be the value of T„ at rj„j. This was completely arbitrary and more 
complicated procedures could be performed to get better fits. Though the plots are 
not identical, they clearly have the same trends, and agreement is within measurement 
error. Using the T3 as a base function this experimental data indicates that T„(r) ^ 
^T^{r). Or, in terms of the unnormalized moments 5'„(r) ^ |^'S'3(S'2)''"~^''/^. The 
coefficients |^ behave in a different manner from their even counterparts, the a„. The 
two sets of coefficients are plotted in Fig. 5.6 along with a dotted line representing 
Gaussian values. Where the a„ are following the Gaussian prediction quite nicely, the 
hn/h^ behave as an almost perfect exponential. 

The experimental results lead to the conclusion that the higher order moments 
can be expressed approximately in terms of the two lowest order moments as 



Sm{r) 



an{S2{r)) 2 : n even. 



dmSs{r){S2{r))'^ 



m odd. 



(5.4) 
(5.5) 
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Figure 5.6: The multiplicative constants a„ and &n/&3 for the data set (c) in Table 
4.1. The dotted line corresponds to the exact values of a purely Gaussian distribution 
given by Eq. 5.3. 

where dm = bm/b^. We are now in a position to draw some conclusions about scale 
invariance in the inertial range. In the inertial range it is clear that 5*3 (r) oc r as 
the exact result from chapter 4 shows. Further, the range 5*2 (r) is not expected to 
deviate in an inertial range from it's r'^^^ behavior. Assuming these two scaling laws 
hold in an inertial range and assuming also that the above results can be extended 
into the inertial range yields 

Sn{r) oc r3 : Vn. (5.6) 

This is precisely the statement of the K41 theory for r in an inertial range. Assuming, 
then, that the extensions and assumptions made are correct, the inertial range of 2D 
turbulence should behave in a scale invariant manner. 

One might ask about the universality of the coefficients a„ and dn- Recall that one 
of the predictions of K41 is that these dimensionless numbers should be independent 
of any external parameters, such as a or F which govern the turbulence. To test this 
the procedure presented above is performed on data set (d) of Table 4.1 as well as a 
set of data taken using a square magnet array instead of a Kolmogorov array. Similar 
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to the earlier data sets, the even moments display strong Gaussian characteristics. 
The underlying structure of the odd normalized moments has changed, but they are 
still marginally collapsible by extracting an arbitrary constant. The odd T„/6„ are 
shown in Fig. 5.7(a) and (b) for the strongly damped Kolmogorov flow (case (d) in 
Table 4.1) and square array respectively. In Fig. 5.8 the coefficients for all of the data 
sets are displayed simultaneously. Note that the measured values do not significantly 
differ from one data set to the next, the a„ remain close to the Gaussian values given 
by Eq. 5.3, and the dn are exponential. This indicates that the coefficients in the 
K41 theory are indeed universal as predicted. 
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Figure 5.7: Tn{r)/hn for odd n for (a) case (d) in Table 4.1 and (b) a run of the e-m 
cell with a square array. 
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Figure 5.8: The measured a„ and dn{= bn/b^) for three data sets using different 
and different types of forcing. 



Chapter 6 
Conclusion 



The inverse energy cascade of 2D turbulence as it occurs in the e-m cell has been 
measured and quantified in the preceding chapters. Clearly an inverse cascade exists 
and it's energy distribution and extent behave as predicted by dimensional analysis. 
The energy flow, when homogenous enough to be accurately compared with theory 
agrees almost perfectly with exact predictions made using the 2D incompressible 
Navier-Stokes equation. However, at no time is this energy flow inertial. 

That the inertial range and inverse cascade range are not coincident for a lin- 
early damped 2D fluid is perhaps the most important result in the thesis, from an 
experimental and numerical stand point. This is because previous 2D turbulence ex- 
periments and simulations did not attempt to differentiate between the two ranges, 
and merely assumed that a —5/3 range in E{k) must be accompanied by the nec- 
essary linear range in 5*3 (r) (or Sj^^ depending on if the turbulence was isotropic or 
not). This is now known not to be the case and casts doubt on some of the earlier ex- 
perimental results. Finally, the inverse cascade range measured in the e-m cell is not 
scale invariant, though there is some evidence that it might become so if an inertial 
range ever developed. 

This is perhaps as far as experimental results on the inverse cascade can be taken 
in the current incarnation of the e-m cell. The lack of an inertial range being the most 
signiflcant limiting feature. To get an idea of what would be needed to eliminate this 
feature we can use the knowledge that has been gained that the inertial range grows 
as the geometric average of the injection and outer scale (Eq. 4.14). If a decade of 
inertial range is desired for accurate measurements of scaling behavior then by Eq. 
4.14 the system size would have to be two decades larger than the injection scale, 
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that is rout = 100 x r,„j . To get an idea what this would mean in the e-m cell, recall 
that the Kolmogorov magnets produced an injection scale of 0.6 cm. Thus the e-m 
cell frame would need to be 60 cm across to support a large enough system to have 
a decade of inertial range while maintaining homogeneity. Moreover, the outer scale 
would have to be forced higher by either increasing einj, which is risky because velocity 
fluctuations would grow and sacrifice incompressibility, or by reducing a. Since a has 
almost reached it's asymptotic limit for the weakly damped cases considered in this 
thesis, a partial vacuum must be employed to achieve lower a. The author need not 
emphasize the difficulty in creating a half meter sized soap bubble, balancing it with 
respect to gravity, and finally putting the whole system in a partial vacuum. 

Some sort of happy medium might be reached by marginally increasing the system 
size, say by a factor of two, and reducing the magnet size a bit. Some inertial 
range would exist, though hardly enough to draw conclusions about scaling behavior. 
However, if use of the e-m cell in inverse cascade investigations is to continue such 
compromises must be made to obtain an inertial range. 



Appendix A 

Particle Tracking Velocimetry: 
Program Listing 

Included in this appendix is the program code, piv_chk.cpp, for the particle tracking 
routine presented in chapter 2. It has been written in the C program language for 
no other reason than personal preference. The program has been successfully com- 
piled with both the GNU C compiler and the Microsoft compiler. No performance 
enhancement was found through use of different compilers. Use of the code after 
compilations is done from the command line by 

piv_chk fool.tif foo2.tif outfoo.vec 

where "fool.tif and "foo2.tif' are the first and second tif images to be compared 
respectively and "outfoo.vec" is the output file holding the particle positions and dis- 
placements. In "outfoo.vec" the origin of the coordinate system is assumed to be at 
the upper left corner of the tif image, with the x direction denoting the vertical coor- 
dinate in the image and the y direction denoting the horizontal coordinate. Values of 
X increase as one goes down the tif image and values of y increase as one moves right. 
If a matched particle set is found at position (xi, yi) in the first image and (x2, ^2) in 
the second image it is saved in "outfoo.vec" as: 

yt "^yi 

where = ^/2{xi + X2), C,y = l/2(yi + 2/2), Ux = X2 — xi and Uy = y2 — yi- All units 
in the output file are in pixels. 
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The various parameters in the routine that determine the characteristics of the 
particles, the search radius, the correlation box size and other parameters are set by 
^define statements at the beginning of the routine. These are commented in the 
source code for clarity. 

piv_chk.cpp: 

#iiiclude <stdio.h> 

#iiiclude <stdlib.h> 

#iiiclude <conio.h> 

#include <math.h> 

/*****Control Paraineters*****/ 

/*The three #define statements below determine the size and properties 
of the pictures. For example, if your picture is (640x480) with 30 
header bits prefacing the contiguous data block in the .tif file, then 
header is 30, row is 480, and col is 640. (use 8 bit .tif)*/ 

#define col 768 
#define row 480 
#define header 234 

/*The two variables below set your background subtraction properties, 
abox is the size of the box to take a local average over and must 
be odd. thresh is the multiple of the background to be subtracted 
from the picture for the purposes of particle identification. Once 
the background has been subtracted, any contiguous group of points 
in the picture with pixel values greater than zero is a candidate 
to be a particle.*/ 

#define abox 21 
#define thresh 1.1 

/*minsize (maxsze) is the minimum (maximum) rms of a candidate group 
of points intensity distribution for that group to be labeled a 
particle. One can think of this as the particle size in pixels.*/ 



#define minsize .25 
#define maxsze 600 
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/*The three variables below are the meat of the routine. 

srad is the distance in pixels to search for the particle from 

frame one to frame two. Set this as low as is reasonable. 

If you know the particles travel no farther than 5 pixels from frame 

one to frame two set srad to 5. cbox is the correlation box size. 

This box must be big enough to contain at least 4 neighbors of 

a particle and must be odd. cthresh is the lower bound of the 

correlation. Any correlation above cthresh will be considered for 

a match, and saved.*/ 

#define srad 15 
#define cbox 17 
#define cthresh 0.5 

/*Most bad interogations happen near the boundaries. 
Set brdr to approx cbox/2 to eliminate these.*/ 

#define brdr 12 

/*The below just allocates memory, maxnum is the maximum # of 
particles in a picture and maxsize is the maximum number of 
contiguous pixels in the candidate block. These are set 
unreasonably high since computer memory is cheap.*/ 

#define maxnum 30000 
#define maxsize 10000 

/*Do not go below here unless you have strongly correlated particle 
motion. If you do not have strongly correlated motion cvstat should 
be zero, which turns the below variables off.*/ 

#define cvstat 
#define nbrstat 2 
#define nrad 15 
#define cvthresh 1.25 
#define Ibnd .6 

/*****Declare Global Variables*****/ 
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struct connection{ 
float val; 
int pic; 

}; 

struct object{ 
float size; 
float x,y; 
float theta; 
int numnbr ,numcan; 
struct connection *nbr; 
struct connection *can; 
int status; 

>; 

/*****Declare Functions*****/ 

void sort (struct connection *,int); 

int mutmax(struct object *, struct object *, int, int); 
void background (unsigned char **, float **,int); 

struct object findparts (unsigned char **, unsigned char **,int,int); 

void connect (struct object *, struct object * , int , int , int) ; 

float correl (unsigned char **, unsigned char **, int); 

void getbox (unsigned char **, unsigned char **, int , int , int) ; 

void clean(struct object *, struct object *, int, int); 

int check_vec (struct object *, struct object *, int, int); 

/*****r/[ain Function*****/ 

void inain(int argv, char *argc[]) 
{ 

FILE *finl, *fin2, *fout; 
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int i,j,k,l,m,n; 
int xO,yO,r,s,t; 
int nl,n2,flag; 
float val; 
float x,y,u,v; 

float thetaO , thetal , thetadif f ; 
float vorticity; 
float distl ,dist2 ,xdif f ,ydif f ; 
double rem.pos; 

unsigned char **picl,**pic2,**bfrl,**bfr2,*hdr; 
float **inean; 

struct object *listl, *list2; 
if (argv<4){ 

printf ("\nsyntax: piv <tif file #1> <tif file #2> <vec file>"); 
exit(O) ; 

> 

/*****Dpen tif and output files*****/ 

if ( (f inl=f openCargc [1] , "rb") )==NULL) { 
printf ("Could not open y„s" , argc [1] ) ; 
exit(O) ; 

} 

if ( (fin2=f openCargc [2] , "rb" ) ) ==NULL) { 
printf ("Could not open "/.s" , argc [2] ) ; 
exit(O) ; 

> 

if ( (f out=f open(argc [3] , " w" ) ) ==NULL) { 
printf ("Could not open "/oS" , argc [3] ) ; 
exit(O) ; 

} 
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/*****Allocate some memory*****/ 
hdr = new (unsigned char [header] ) ; 
picl = new(unsigned char * [row] ) ; 
pic2 = new(unsigned char * [row] ) ; 
bfrl = new(unsigned char * [row] ) ; 
bfr2 = new (unsigned char * [row] ) ; 
mean = new(float * [row] ) ; 

f or (i=0 ; i<row; i++) { 

picl[i] = new(unsigned char [col]); 
pic2[i] = new(unsigned char [col]); 
bf rl [i] = new(unsigned char [col]); 
bfr2[i] = new(unsigned char [col]); 
mean[i] = new(f loat [col] ) ; 

> 

/*****Read picture files*****/ 
printf ("Reading files\n"); 

if (fread(hdr, sizeof (unsigned char) .header ,finl) !=header)-C 
printf ("Error stripping header from "/oS" , argc [1] ) ; 
exit(O) ; 

> 

if (fread(hdr, sizeof (unsigned char) .header ,fin2) !=header){ 
printf ("Error stripping header from 7oS" , argc [2] ) ; 
exit(O) ; 

} 

f or (i=0 ; i<row; i++) { 

if (f read (picl [i] , sizeof (unsigned char) ,col,finl) !=col){ 
printf ("Error reading °/oS" ,argc [1] ) ; 
exit(O) ; 

} 

if (f read(pic2 [i] , sizeof (unsigned char) , col ,f in2) ! =col) { 
printf ("Error reading °/oS" ,argc [2] ) ; 
exit(O) ; 

> 
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/*****Begin background subtraction*****/ 
printf ("\nSubtracting background") ; 
backgroundCpicl ,mean, abox) ; 
f or (i=0 ; i<row ; i++) { 
for(j=0; j<col; j++){ 

if (thresh*mean[i] [j] >= picl[i][j]) bfrl[i][j] = 0; 
else{ 

val = mean[i] [j] ; 
rem = modf (val ,&pos) ; 
if (rem > .5) pos++; 
val =(float)pos; 

bfrl[i][j] = picl[i][j] - (unsigned char)val; 

} 

} 

} 

background (pi c2 , mean, abox) ; 
f or (i=0 ; i<row; i++) { 
for(j=0; j<col; j++){ 

if (thresh*mean[i] [j] >= pic2[i][j]) bfr2[i][j] = 0; 
else{ 

val = mean[i] [j] ; 
rem = modf (val ,&pos) ; 
if (rem > .5) pos++; 
val =(float)pos; 

bfr2[i][j] = pic2[i][j] - (unsigned char)val; 

} 

} 

} 

little housekeeping*****/ 
delete [] mean; 

listl = new(struct object [maxnum] ) ; 
list2 = new(struct object [maxnum] ) ; 
nl = 0; 
n2 = 0; 
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/*****Begin Finding particles*****/ 
printf ("\nFinding particles"); 
f or(i=0; i<row; i++){ 
for(j=0; j<col; j++){ 
if (bfrlEi] [j] > 0){ 

listl [nl] =f indparts (bf rl ,picl , i , j ) ; 

if (listl [nl] . size < maxsze && listl [nl] . size > minsize/2 ) nl++ 

} 

if (bfr2[i] [j] > 0){ 

list 2 [n2] =f indparts (bf r2,pic2 , i , j ) ; 

if (list2 [n2] . size < maxsze && list2 [n2] . size > minsize/2 ) n2++ 

} 

} 

> 

little housekeeping*****/ 

delete [] bfrl; 
delete [] bfr2; 

/*****Start finding neighbors and candidates*****/ 

printf ("\n°/od\t°/od\nEstablishing connections" ,nl,n2) ; 
connect (listl , list2 ,nl ,n2 , srad) ; 
bfrl = new(unsigned char * [cbox] ) ; 
bfr2 = new(unsigned char * [cbox] ) ; 
f or(i=0; Kcbox; i++){ 

bf rl [i] = new(unsigned char [cbox]); 

bfr2[i] = new(unsigned char [cbox]); 

for(j=0; j<cbox; j++){ 
bfrl[i] [j] = 0; 
bfr2[i] [j] =0; 

} 
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/*****Calculate correlations for the candidates*****/ 
printf ("\nBeginning correlations") ; 
for(i=0;i<nl;i++){ 

ifCiy.lOO == 0) printf ("."); 

m = listl [i] .numcan; 

if (m==0) continue; 

ifCkbhitO) break; 

val = listl [i] .x; 
rem = modf (val,&pos) ; 
if (rem > .5) pos++; 
xO =(int)pos; 
val = listl [i] .y ; 
rem = modf (val ,&pos) ; 
if (rem > .5) pos++; 
yO =(int)pos; 

getbox(picl ,bf rl ,xO ,yO , cbox) ; 

for(j=0; j<m; j++){ 

n=listl[i] .can[j] .pic; 
val = list2[n] .x; 
rem = modf (val ,&pos) ; 
if (rem > .5) pos++; 
xO = (int)pos; 
val = list2 [n] .y ; 
rem = modf (val ,&pos) ; 
if (rem > .5) pos++; 
yO = (int)pos; 

getbox (pic2 , bf r2 , xO , yO , cbox) ; 

val = correl(bfrl,bfr2,cbox) ; 
listl [i] . can [j] . val = val; 
xO = 0; 

while (list2 [n] . can [xO] .pic != i) xO++; 
list2[n] .can[xO] .val = val; 

} 

> 
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/*****Sort by size of correlations*****/ 

printf ("\nsorting") ; 
for(i=0;i<nl;i++){ 

m = listl [i] .numcan; 

sort (listl [i] . can,m) ; 

m = listl [i] .numnbr; 

sortdistl [i] .nbr,in) ; 

} 

for(i=0;i<n2;i++){ 
m = list2 [i] .numcan; 
sort (list2 [i] . can,ni) ; 
m = list2 [i] .numnbr ; 
sort (list2 [i] .nbr,m) ; 

> 

/*****Start matching*****/ 
printf ("\nmatching") ; 
delete [] picl; 
delete [] pic2; 

f or ( val= . 9 ; val>cthresh ; val-=0 . 01) { 
for(i=0;i<nl;i++){ 

if (listl [i] .numcan==0) continue; 
if (listl [i] . status==l) continue; 
m = listl [i] .numcan - 1; 
n = listl [i] . can [m] .pic; 

if (mutmaxdistl , list2, i ,n)==l && listl [i] . can [m] .val>=val){ 
if(cvstat == 1 && val < lbnd){ 

if (check_vec (listl , list2 , i >n)==l){ 
listl [i] . status = 1; 
list2 [n] . status = 1; 
cleandistl , lists , i ,n) ; 

> 
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else if (check_vec (listl , list2, i ,n)==0){ 
listl[i] .numcan — ; 
list2[n] .numcan — ; 

> 

} 

else{ 

listl [i] . status = 1; 
list2 [n] . status = 1 ; 
clean(listl,list2,i,n) ; 

} 

} 

} 

} 

n=0; 

/*****Print results*****/ 
for(i=0;i<nl;i++){ 

if (listl [i] .status == 1){ 

j = listl [i] .numcan-l ; 
m = listl [i] . can[j] .pic; 

if (listl [i] . size < minsize || list2 [m] . size < minsize) continue; 
X = (list2[ni].x + listl [i] .x)/2; 
y = (list2[m].y + listl [i] . y) /2 ; 
u = list2[iii].x - listl[i].x; 
V = list2[iii].y - listl[i].y; 

if (x > brdr && x <= row-brdr && y > brdr && y<= col-brdr){ 
f printf (font , ""/.f \ty.f \ty.f Xt/^f \n" , y , x , v ,u) ; 
n++; 

} 

} 

} 

printf ("\n%d" ,n) ; 
/*****Soine final housekeeping*****/ 
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delete [] listl; 
delete [] list2; 

} 

/*This function attempts to loosen the correlation limit by comparing 
candidate motions with previously matched neighbors.*/ 

int check_vec (struct object *listl , struct object *list2,int nl, int n2) 

■c 

int numl,num2; 
int nbrs , matched; 
int i,j,m,n; 
float x,y,u,v,val; 
float xm,ym,um,vm; 
float xp,yp,up,vp; 
float jtr; 

matched=0 ; 

nbrs = listl [nl] .numnbr; 
xp = (Iist2[n2] .x+listl [nl] 
yp = (list2 [n2] .y+listl [nl] 
up = (Iist2[n2] .x-listl [nl] 
vp = (Iist2[n2] .y-listl [nl] 

f or (i=0 ; Knbrs ; i++) { 
m=listl[nl] .nbr[i] .pic; 
if (listl [m] . status == 1){ 

j=listl[m] .numcan-1; 

n=listl[m] .can[j] .pic; 

X = (list2[n] .x+listl [m] .x)/2 - xp; 

y = (list2[n] .y+listl [m] .y)/2 - yp; 

val = sqrt(x*x + y*y) ; 
if (val < nrad && x+xp > brdr && x+xp<=row-brdr 

&& y+yp > brdr && y+yp <= col-brdr) matched++; 

} 

} 



.x)/2; 
.y)/2; 

• x); 

• y); 



if (matched < nbrstat) return(2) ; 
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xm=0 ; ym=0 ; um=0 ; vm=0 ; 
f or (i=0 ; i<nbrs ; i++) { 
m=listl [nl] .nbr[i] .pic; 

if (listl [m] . status == 1){ 
j=listl[m] .numcan-1; 
n=listl[m] .can[j] .pic; 
X = (list2[n] .x+listl [m] .x)/2 - xp; 
y = (list2[n] .y+listl [m] .y)/2 - yp; 
val = sqrt(x*x + y*y) ; 

if (val < nrad && x+xp > brdr && x+xp<=row-brdr 

y+yp > brdr && y+yp <= col-brdr){ 
xm += (list2[n] .x+listl [m] .x)/(2*(float)matched) ; 
ym += (list2 [n] .y+listl [m] .y)/(2*(f loat)matched) ; 
urn += (list2[n] .x-listl [m] .x)/(float)niatched; 
vm += (list2 [n] .y-listl [m] .y)/(f loat)niatched; 

} 

> 

> 

jtr=0; 

f or (i=0 ; Knbrs ; i++) { 
m=listl [nl] .nbr[i] .pic; 

if (listl [m] . status == 1){ 
j=listl[m] .numcan-1; 
n=listl[m] .can[j] .pic; 
X = (list2[n] .x+listl [m] .x)/2 - xp; 
y = (list2[n] .y+listl [m] .y)/2 - yp; 
val = sqrt(x*x + y*y) ; 

if (val < nrad && x+xp > brdr && x+xp<=row-brdr 

&& y+yp > brdr && y+yp <= col-brdr){ 
u = (list2 [n] . x-listl [m] . x) - urn; 
V = (list2 [n] . y-listl [m] . y) - vm; 
jtr += (u*u+v*v) / (f loat)matched; 

} 
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> 

} 



jtr=sqrt ( jtr) ; 
u = up-um; 
V = vp-vm; 

val = sqrt (u*u + v*v) ; 

if (val < cvthresh*jtr) return(l) ; 

else return(O) ; 



/♦Cleans the objects i and j from any candidate list since they 
have presumably been matched.*/ 

void clean(struct object *listl , struct object *list2,int nl, int n2) 
{ 

int numl,num2; 

int i,j,k,l,m; 

struct connection temp; 

numl = listl [nl] .numcan - 1; 



if(numl > 0){ 

f or (i=0 ; i<numl ; i++) { 

j = listl [nl] . can[i] .pic; 
m = list2 [j] .numcan; 
temp = list2 [j] . can [m-1] ; 
k=0; 

while(list2[j] .can[k] .pic != nl) k++; 
list2 [j] . can [k] = temp; 
list2[j] .numcan — ; 
m— ; 

sort (list2 [j] .can,m) ; 

> 

} 



num2 = list2 [n2] .numcan - 1; 
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if(iiuiii2 > o){ 

f or(i=0; i<num2; i++){ 

j = list2 [n2] . caii[i] .pic ; 
m = listl [j] .numcan; 
temp = listl [j] . can [m-1] ; 
k=0; 

while (list 1 [j] . can [k] .pic != n2) k++; 
listl [j] . can [k] = temp; 
listl [j] .numcan — ; 
m— ; 

sort (listl [j] . can,m) ; 

} 

> 

} 

/*Returns a 1 if the maximum correlation of listl [nl] and Iist2[n2] 
point to one another. otherwise. This routine assumes we 
have already sorted the connections using "sort".*/ 

int mutmax(struct object *listl , struct object *list2,int nl.int n2) 
{ 

int numl,num2; 
int i , j ; 

numl = listl [nl] .numcan - 1; 
num2 = list2 [n2] .numcan - 1; 
i = listl [nl] . can [numl] .pic; 
j = list2 [n2] . can[num2] .pic; 
if(i == n2 && j == nl) return ( 1) ; 
else return(O) ; 

} 

/*Sorting routine used in object lists. This is used to sort the 
nbr connections in order from closest to farthest away and the can 
connection from lowest correlation to highest.*/ 
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void sort (struct connection *ptr,int num) 
{ 

int i , j ; 

struct connection temp; 
for(j=l; j<num; j++){ 

temp = ptr[j] ; 

i=j-l; 

while(i>=0 && ptr[i].val > temp.val){ 
ptr[i+l]=ptr[i] ; 
i~; 

} 

ptr[i+l] = temp; 

} 

} 

/*Finds all the parts of an particle given that there is a bright spot 
at xO,yO. Return the value of the centroid and the rms.*/ 

struct object findparts (unsigned char **ptrl .unsigned char **ptr2, 

int xO,int yO) 

{ 

int i,m=l,n=l; 
int *x,*y; 
int j ,k,val; 

unsigned char *intensity; 
int brght=0; 
float tx,ty,std; 
struct object out; 

X = new(int [maxsize] ) ; 

y = new (int [maxsize] ) ; 

x[0] =xO; 

y[0] =yO; 

ptrl [xO] [yO]=0; 



while (n<maxsize) { 
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f or(i=0; i<m; i++){ 
if(x[i] - 1 >= 0){ 

if (ptrl[x[i]-l] [y[i]] > 0){ 
x[n] = x[i] - 1; 
y[n] = y[i] ; 
ptrl[x[n]] [y[n]]=0; 
n++; 

} 

} 

if(x[i] + 1 < row){ 

if (ptrl[x[i]+l] [y[i]] > 0){ 
x[n] = x[i] + 1; 
y[n] = y[i] ; 
ptrl[x[n]] [y[n]]=0; 
n++; 

} 

} 

if(y[i] + 1 < col){ 

if (ptrl[x[i]] [y[i]+l] > 0){ 
x[n] = x[i] ; 
y[n] = y[i]+l; 
ptrl[x[n]] [y[n]]=0; 
n++; 

> 

} 

if(y[i] - 1 >= 0){ 

if (ptrl[x[i]] [y[i]-l] > 0){ 
x[n] = x[i] ; 
y[n] = y[i]-l; 
ptrl[x[n]] [y[n]]=0; 
I1++; 

} 

} 

> 
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if (n==m) break; 
else m=n; 

} 

intensity = new(unsigned char[n]); 

for(i=0;i<n;i++){ 

intensity [i] = ptr2[x[i]] [y[i]] ; 
brght += (int) intensity [i] ; 

} 

out . x=0 ; 
out . y=0 ; 

f or(i=0; i<n; i++){ 

out.x += (float) intensity [i] * (float)x [i] / (float)brght ; 
out.y += (float) intensity [i] * (float)y [i] / (float)brght ; 

> 

std = 0; 

for(i=0;i<n;i++){ 

tx = (f loat)x [i] -out .x; 
ty = (f loat)y [i] -out .y; 

std += (float) intensity [i] * (tx*tx+ty*ty) / (float)brght ; 

} 

std = sqrt (std) ; 
out. size = std; 
out . numnbr = ; 
out.numcan = 0; 
out. status = 0; 

delete [] intensity; 
delete [] x; 
delete [] y; 
return (out) ; 
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/♦Creates float ** mean which contains the value of the mean 
of a box of sizexsize around each point in pic.*/ 

void background (unsigned char **pic, float **back,int size) 
{ 

int i , j ; 

int x,y,val=0; 

int half = (size-l)/2; 

float mean; 

mean=0 ; 

f or (i=0 ; i<=half ; i++) { 
for(j=0; j<=half ; j++){ 

mean = (val*mean +(f loat)pic [i] [j] )/( (float) val+1) ; 
val++; 

} 

} 

back[0] [0]=mean; 

for(j=0; j<row; j+=2){ 

if(jy.64==0) printf ("."); 
for(i=l;i<col;i++){ 
y=i-half-l ; 

for (x= j -half ; x<= j +half ; x++) { 

if (y < II X < II x >= row) continue; 

mean = (mean*val - (f loat)pic [x] [y] )/( (float) val-1) ; 

val — ; 

} 

y=i+half ; 

for (x= j -half ; x<= j +half ; x++) { 

if (y >= col II X < II X >= row) continue; 

mean = (mean*val + (f loat)pic [x] [y] )/( (float) val+1) ; 

val++; 

} 
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back[j] [i]=mean; 



X = j-half; 

f or (y=col-l-half ; y<col ; y++) { 
if (x < 0) continue; 

mean =(mean*val - (f loat)pic [x] [y] ) / ( (f loat)val-l) ; 
val — ; 



X = j+l+half; 

f or (y=col-l-half ; y<col ; y++) { 
if (x >= row) continue; 

mean =(mean*val + (f loat)pic [x] [y] )/( (float) val+1) ; 
val++; 

> 

back[j+l] [col-1] = mean; 
f or(i=col-2; i>=0; i — ){ 
y=i+half+l; 

f or (x=j +l-half;x<= j+l+half ;x++){ 

if (y >= col II X < II X >= row) continue; 

mean = (mean*val - (f loat)pic [x] [y] )/( (float) val-1) ; 

val — ; 

} 

y=i-half ; 

for (x= j + l-half ; x<= j + l+half ; x++) { 

if (y < II X < II X >= row) continue; 

mean = (mean*val + (f loat)pic [x] [y] )/( (float) val+1) ; 

val++; 

} 



back[j + l] [i]=mean; 

> 
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X = j+l-half; 
f or(y=0;y<=half ;y++){ 
if (x < 0) continue; 

mean =(mean*val - (f loat)pic [x] [y] )/( (float) val-1) ; 
val — ; 

> 

X = j+2+half; 
f or(y=0;y<=half ;y++){ 
if (x >= row) continue; 

mean =(mean*val + (f loat)pic [x] [y] )/( (float) val+1) ; 
val++; 

> 

if((j+2) >= row) continue; 
back [j +2] [0] = mean; 

} 

} 

/*Finds all the neighbors and candidates for a particle and then 
stores this info in the appropriate spot in the lists.*/ 

void connect (struct object *listl , struct object *list2, 

int nl,int n2,int maxdist) 

{ 

int xO,yO; 

int i , j ; 

int x,y; 

int ml,m2; 

int **picl , **pic2 ; 

int *templ,*temp2; 

float val ,dist ,xdif f ,ydif f ; 

double rem,pos; 



picl 
pic2 



= new (int * [row] ) ; 
= new (int * [row] ) ; 
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f or (i=0 ; i<row; i++) { 

picl[i] = new(int [col] ) ; 
pic2[i] = new(int [col] ) ; 
for(j=0; j<col; j++){ 

picl[i] [j] = -1; 

pic2[i] [j] = -1; 

} 



for(i=0;i<nl;i++){ 
val = listl [i] .x; 
rem = modf (val,&pos) ; 
if (rem >.5) pos++; 
X = (int)pos; 
val = listl [i] .y ; 
rem = modf (val, &pos) ; 
if (rem >.5) pos++; 
y = (int)pos; 
picl [x] [y] = i; 



for(i=0;i<n2;i++){ 
val = list2 [i] . x; 
rem = modf (val, &pos) ; 
if (rem >.5) pos++; 
X = (int)pos; 
val = list2 [i] .y; 
rem = modf (val ,&pos) ; 
if (rem >.5) pos++; 
y = (int)pos; 
pic2[x] [y] = i; 



tempi = newdnt [maxdist*maxdist] ) ; 
temp2 = new(int [maxdist*maxdist] ) ; 
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f or (xO=0 ; xO<row ; xO++) { 

if(xO 7, 64 == 0) printf (" . ") ; 
for (yO=0 ; yO<col ; yO++) { 

if (picl [xO] [yO] == -1 && pic2 [xO] [yO] ==-1) continue; 
if (piclExO] [yO] != -1){ 
inl=0; 
in2=0; 

f or (i=-maxdist ; i<=maxdist ; i++) { 
X = xO + i ; 

if (x < II X >= row) continue; 
f or ( j =-niaxdist ; j <=niaxdist ; j ++) { 

y = yo + j; 

if (y < II y >= col) continue; 

if((i*i + j*j)> maxdist*niaxdist) continue; 

if (pic2[x] [y] != -1){ 

temp2 [m2] = pic2 [x] [y] ; 

m2++; 

} 

if (picl [x] [y] != -1){ 

if (x==xO && y==yO) continue; 
tempi [ml] = picl[x][y]; 
ml++; 

> 

> 

} 

listl [picl [xO] [yO] ] .nbr = new(struct connection [ml] ) ; 

for(i=0;i<ml;i++){ 

xdiff = listl [tempi [i]] .X - listl [picl [xO] [yO] ] .x; 
ydiff = listl [tempi [i]] .y - listl [picl [xO] [yO]] .y; 
dist = sqrt(xdiff *xdiff+ydiff *ydiff ) ; 
listl [picl [xO] [yO]] .nbr [i] .pic = tempi [i] ; 
listl [picl [xO] [yO]] .nbr[i] .val = dist; 
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listl [picl [xO] [yO]] .numnbr = ml; 

listl [picl [xO] [yO] ] . can = new(struct connection [m2] ) ; 

f or(i=0; i<m2; i++){ 

listl [picl [xO] [yO] ]. can [i] .pic = temp2 [i] ; 

} 

listl [picl [xO] [yO]] .numcan = m2; 



if (pic2[x0] [yO] != -1){ 
inl=0; 
m2=0; 

f or (i=-maxdist ; i<=maxdist ; i++) { 
X = xO + i ; 

if (x < II X >= row) continue; 
f or ( j =-maxdi St ; j <=maxdi st ; j ++) { 

y = yo + j; 

if (y < II y >= col) continue; 

if((i*i + j*j)> maxdist*niaxdist) continue; 

if (picl [x] [y] != -1){ 

temp2[m2] = picl[x] [y] ; 

m2++; 

> 

if (pic2[x] [y] != -1){ 

if (x==xO && y==yO) continue; 
tempi [ml] = pic2[x][y]; 
ml++; 

} 

} 

} 

list2 [pic2 [xO] [yO] ] .nbr = new(struct connection [ml] ) ; 
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for(i=0;i<ml;i++){ 

xdiff = list2 [tempi [i]] .X - list2 [pic2 [xO] [yO] ] .x; 
ydiff = list2 [tempi [i]] .y - list2 [pic2 [xO] [yO] ] .y; 
dist = sqrt (xdiff *xdiff + ydiff *ydiff) ; 
Iist2[pic2[x0] [yO]] .nbr [i] .pic = tempi [i] ; 
Iist2[pic2[x0] [yO]] .nbr[i] .val = dist; 

} 

list2 [pic2 [xO] [yO]] .numnbr = ml; 

list2 [pic2 [xO] [yO] ] . can = new(struct connection [m2] ) ; 

f or(i=0; i<m2; i++){ 

Iist2[pic2[x0] [yO]] .can[i] .pic = temp2[i] ; 

} 

list2 [pic2 [xO] [yO]] .numcan = m2; 

} 

> 

} 



delete [] picl; 

delete [] pic2; 

delete [] tempi; 

delete [] temp2; 



/*Returns the correlation number between two arrays of size x size.*/ 

float correl (unsigned char **ptrl , unsigned char **ptr2,int size) 
{ 

int i , j ; 

float meanl,mean2; 
float stdl,std2,cor; 
meanl=0 ; 
mean2=0; 
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f or (i=0 ; i<size ; i++) { 
for(j=0; j<size; j++){ 

meanl += (float)ptrl [i] [j] / (float) (size*size) ; 
inean2 += (float)ptr2 [i] [j] / (float) (size*size) ; 

} 

} 

f or (i=0 ; i<size ; i++) { 
for(j=0; j<size; j++){ 

stdl += (((float)ptrl [i] [j] - meanl) *( (float)ptrl [i] [j] - meanl)) 

/(float) (size*size) ; 
std2 += (((float)ptr2[i] [j] - mean2) * ( (f loat)ptr2 [i] [j] - mean2)) 

/(float) (size*size) ; 
cor += (((float)ptrl [i] [j] - meanl) *( (float)ptr2 [i] [j] - mean2)) 
/ (float) (size*size) ; 

} 

} 

cor /=(sqrt(stdl)*sqrt(std2)) ; 
return (cor) ; 

} 

/*Gets a box from picl centered at xO,yO and stores it in ptrl. 
ptrl should be at least size x size and size must be odd! ! ! .*/ 

void getbox (unsigned char **picl, unsigned char **ptrl, 

int xO,int yO,int size) 

{ 

int i,j,x,y; 

int half = (size-l)/2; 

f or (i=0 ; i<size ; i++) { 
X = xO + i - half ; 
if (x<0 I I x>=row){ 
for(j=0; j<size; j++){ 
ptrlEi] [j] = 0; 

} 

continue ; 

> 
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f or ( j=0 ; j<size ; { 
y = yO + j - half; 
if (y<0 I I y>=col){ 
ptrlEi] [j] = 0; 
continue ; 

} 

ptrlEi] [j] = picl[x] [y] ; 

> 
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